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ABSTRACT 
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CHAPTER  I 


BACKGROUND  AND  RESEARCH  OBJECTIVE 


1.1  Introduction 

Both  partial  differential  equations  and  topology  are  closely  related 
to  classical  mechanics,  the  former  through  Hamilton-Jacobi  theory  and 
the  latter  through  global  analysis.  In  Hamilton-Jacobi  theory,  a 
configuration  space  differential  equation  is  given  a  phase  space  represen¬ 
tation.  In  global  analysis,  solutions  of  differential  equations  are 
studied  geometrically  in  phase  space.  In  this  investigation, the  asymptotic 
series  solution  of  a  particular  partial  differential  equation,  the 
reduced  Helmholtz  wave  equation,  is  determined  using  topological 
considerations  in  an  essential  way. 


1.2  Background 

Wave  propagation  at  a  definite  frequency,  w,  is  commonly  represented 
by  a  partial  differential  equation  of  the  form 


„2 ,  —  . ,  f  (r,w)  81 

V  i{i(r,t)  =  s  l  L  — , 

c  3t 


(1.2.1) 


where  ijj(r,t)  is  the  wave  function,  r  refers  to  the  spatial  coordinates, 
t  is  the  time,  f(r,uj)  is  called  the  profile  characterizing  the 
inhomogeneity  of  the  medium  and  c  is  the  phase  velocity  when  the 
medium  is  homogeneous,  i.e.,  f(r,w)  =1.  If  the  medium  is  non-dispersive, 
the  profile  is  frequency  independent,  i.e.,  f (r ,to)  = f (r) .  For  waves 
transmitted  at  a  particular  frequency,  the  time  dependence  may  be 
represented  by  exp{iwt).  If,  for  a  non-dispersive  medium. 


if>(r,t)  ■  iHr)exp{iwt} 
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Equation  (1.2.1)  becomes 

V2ip(r)  +tf (r)ip(r)  =  0  ,  (1.2.2) 

where 

t  =w2/c2  .  (1.2.3) 

This  is  the  reduced  Helmholtz  wave  equation,  a  second  order  partial 
differential  equation  with  a  variable  coefficient.  Although  no  general 
technique  exists  for  solving  such  equations  [42],  various  techniques  do 
exist  for  constructing  solutions  valid  under  specific  assumptions  [11,  46]. 
One  such  technique,  valid  at  high  frequencies,  is  the  approximation  of 
iji(r)  by  an  asymptotic  series. 

In  the  asymptotic  series  solution  of  Equation  (1.2.2),  it  is  assumed 

that 

t|j(r)  -  exp{ it4> (r)  }  l  a ,(r)x  k  ,  (1.2.4) 

k=0  K 

as  t  ■*■*>.  Substitution  of  Equation  (1.2.4)  into  Equation  (1.2.2), 
regrouping  by  powers  of  t,  then  equating  the  coefficients  of  the  powers 
of  t  to  zero  yields 

I  (V<f>)|  2  -  f(7)  =  0  (1.2.5) 

and 

2(V<j>)  •  (Vak)  + (V2<J.)ak  =  -V2ak_1;  k- 0,1,2,...;  a^-0.  (1.2.6) 

Equation  (1.2.5)  is  the  eikonal  equation;  Equation  (1.2.6)  is  the 
transport  equation.  ^ (r )  is  obtained  by  solving  the  eikonal  equation 
for  <t>(r),  then  substituting  $(r)  into  the  transport  equation  determines 
the  ak(r)  recursively.  The  first  term  in  the  series,  ao(r)exp{ix<f>(r) ) , 
is  the  JWKB  approximation  [30].  The  remaining  terms  may  be  regarded 
as  higher  order  corrections  to  JWKB  approximation  [52]. 


3 

Two  basic  difficulties  with  this  approach  are  noteworthy.  First 


<j)(r)  must  be  determined  from  a  first  order  nonlinear  partial  differential 
equation.  As  no  general  technique  exists  for  solving  such  equations, 

<j>(r)  itself  may  require  an  approximate  solution.  The  second  difficulty 
involves  singular,  or  turning,  points  and  is  most  clearly  illustrated 
by  considering  the  eikonal  and  transport  equations  in  one  dimension. 


i.e. , 


«&>  -fw-o 


(1.2.7) 


and 


d2a. 


2di  ^!k  +  d!ia  .  .  k=0  i  2 

dx  dx  .  2  ak  ,2  ’ 

dx  dx 


;  a_1  =  o. 


(1.2.8) 


Equation  (1.2.7)  can  be  integrated  to  obtain 
rx 


4>(x)  =  + 


f(s)1/2ds  . 


Substituting  for  <j>(x)  in  Equation  (1.2.8)  determines 
aQ(x)  =  f (x)  1/^2  . 

d<Hxo) 

If  at  x  =  x  ,  — - -  =  0,  then  f(x  )  =  0  and  a  (x  )  is  unbounded.  Points 

o  dx  o  o  o 

such  as  xq  are  called  singular  points.  Points,  x^,  at  which  aQ(x) 

does  not  become  unbounded  are  called  regular  points. 

More  physically,  in  two  or  three  dimensions,  a  finite  source 

produces  a  field  of  waves  which  can  be  described  at  high  frequencies 

by  a  family  of  trajectories.  The  contact  of  neighboring  trajectories 

near  turning  points  leads  to  a  focusing  which  effects  a  theoretically 

infinite  field  intensity  on  an  envelope  defined  by  the  locus  of  turning 

(singular)  points.  Beyond  this  boundary,  in  theory,  no  energy  is 
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propagated.  Heuristically ,  a  caustic  is  the  higher  dimensional  analog 
of  a  turning  point;  analogously,  the  above  asymptotic  expansion  is 
not  valid  on  caustics  [12]. 

The  asymptotic  expansion  technique  has  been  modified  by  several 
investigators  to  extend  its  validity  to  include  caustic  regions.  One 
approach  is  to  regard  the  caustic  as  a  boundary  on  which  a  separate 
asymptotic  expansion  is  determined;  however,  this  boundary  expansion 
must  be  matched  termwise  with  expansions  valid  in  non-caustic  regions  [46] . 
Another  approach  is  to  determine  a  uniformly  valid  asymptotic  series, 
i.e.,  an  asymptotic  series  which  satisfies  the  reduced  Helmholtz  wave 
equation  to  arbitrarily  high  order  in  t  ^  in  regions  which  contain 
caustics.  The  procedure  involves  transforming  a(r)  to  an  appropriate 
special  function,  e.g.,  Bessel,  Hankel,  etc.,  and  determining  the 
asymptotic  series  in  terms  of  the  special  function.  The  technique 
applies  only  to  caustics  of  simple  geometry  [54]. 

An  alternate  formulation  is  an  integral  representation  for  t^(r) , 


i.e.. 


■Kr) 


' 

A(r,p,T)exp{ix<Kr,p)  )dp 

J, 


(1.2.9) 


where 

A(r,p,x)  -  T.  A,  (r ,p)x  k  . 

k=0 

In  such  a  representation  A(r,p,x)  may  be  regarded  as  an  amplitude  and 
4> (r , p)  as  a  phase;  p  is  the  wave  vector,  i.e.,  the  normal  to  the  wave- 
fronts,  the  surfaces  of  constant  phase  [14,  33].  More  physically,  the 
field  at  any  point  is  represented  as  a  superposition  of  incident  plane 
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waves  summed  over  all  directions.  Consequently,  Equation  (1.2.2)  has 
an  asymptotic  solution  of  the  form 


'J'(r)  - 


A(r,p,x)exp{iT<t>(r,p)  }dp  =  0(t”  ) 


(1.2.10) 


as  The  full  asymptotic  series  of  t|>(r)  is  then  the  sum  of  the 

asymptotic  series  of  each  integral 


JJ 


Ak(r,p)exp{iti}i(r,p)}dp 


Maslov  (as  explained  by  Arnold  [10],  Eckmann  and  Seneor  [24]  and 
Guillemin  and  Sternberg  [28])  extended  the  significance  of  this 
representation.  He  noted  that  if  the  wave  vector,  p,  was  regarded  as 
a  momentum,  the  phase,  $(r,p),  could  be  determined  from  considerations 
in  Hamllton-Jacobi  theory,  rather  than  directly  from  the  eikonal  equation. 
Maslov  also  determined  a  formula  for  the  first  term  in  the  asymptotic 
series,  i.e.,  the  asymptotic  limit,  at  regular  points.  The  approach 
of  Maslov  was  generalized  to  include  singular  points  by  Arnold  [2-9]. 

The  entire  asymptotic  series  at  regular  points  of  the  integral 
in  Equation  (1.2.9)  may  be  determined  by  using  the  stationary  phase 
technique  [11].  There  are  two  basic  approaches  to  the  classical 

technique.  One  involves  transforming  the  integral  to  the  form 

' 

exp{iru}dg(u)  , 

. 

where  g(u)  is  the  integral  of  A(r,p,t)  over  the  region  <|>(r,p)  <^u  [32], 

The  asymptotic  behavior  is  then  determined  by  studying  g(u) .  The 
other  approach  is  to  investigate  the  asymptotic  behavior  by  studying 
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the  integral  in  Equation  (1.2.9)  directly  [15].  In  this  research,  the 
latter  approach  is  followed. 

The  algebraic  complexity  of  applying  the  stationary  phase  technique 
to  multiple  integrals  has  been  noted  [19].  Specifically,  explicit 
coordinate  transformations  which  carry  the  integrals  to  tractable  forms 
must  be  determined.  Consequently,  the  technique  is  most  clearly 
illustrated  by  considering  the  simple  integral 
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geometry  by  transforming  $(r,p)  to  an  Airy  function.  The  series  is  then 
the  sum  of  the  Airy  integrals,  which  have  been  tabulated  [1].  In  his 
formulation,  however,  Ludwig  obtained  ifi(r,p)  from  the  eikonal  equation 
in  momentum  space. 

Duistermaat  [22,23],  extending  both  the  work  of  Maslov  and  Ludwig, 
developed  a  formulation  that  superseded  all  previous  approaches.  Maslov 
and  Arnold  [3-10]  realized  that,  at  singular  stationary  points,  <(>(r,p) 
was  topologically  equivalent  to  a  canonical  Thom  potential.  The 
equivalence  is  determined  by  the  geometry  of  the  caustic.  (The  Airy 
function  is  the  lowest  order  Thom  potential.)  Duistermaat,  however, 
realized  that  the  transformation  to  Thom  potentials  led  to  uniformly 
valid  asymptotic  series  for  regions  containing  any  physically  realizable 
caustic  geometry,  thus  extending  Ludwig's  formulation.  (Duistermaat  also 
extended  Maslov's  formula  for  the  asymptotic  limit  to  determine  the  full 
asymptotic  series  at  regular  stationary  points.)  As  in  Ludwig's 
approach,  this  asymptotic  series  is  determined  as  a  sum  of  the  trans¬ 
formed  integrals.  Unlike  Ludwig's  approach,  however,  Duistermaat 's 
algorithm  was  not  entirely  explicit.  Neither  were  the  coordinate 
transformations  carrying  4>(r,p)  to  the  Thom  potentials  explicitly 
given,  nor  were  the  final  canonical  integrals  evaluated. 

1.3  Scope  of  This  Investigation 

In  this  investigation,  phase  space  analysis  is  used  to  determine 
explicitly  the  asymptotic  solution  of  the  reduced  Helmholtz  wave 
equation.  Each  field  point  is  represented  by  the  range,  i.e.,  the 
distance  between  two  points  on  a  horizontal  plane,  and  the  depth.  Hence, 
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two  coordinates  suffice  to  specify  the  spatial  variation  of  ij; (r )  .  For 


computational  purposes,  a  common  technique  is  not  to  represent  the 


medium  as  a  continuum  but  as  a  series  of  horizontal  layers  [17],  Rather 


than  characterizing  the  medium  with  a  single,  depth-dependent  profile. 


each  layer  is  characterized  separately  by  a  linear  approximation  to  the 


profile  valid  only  within  the  layer.  In  this  investigation,  however. 


the  medium  is  characterized  by  a  single,  depth-dependent  profile. 


f(x),  which  is  assumed  to  be  invertible,  i.e.,  x  =  f  (y)  is  defined. 


except  at  a  finite  number  of  isolated  points.  (At  such  points,  the 


algorithm  detailed  in  Chapter  VII  does  not  apply.)  Thus,  the  reduced 


Helmholtz  equation  considered  here  is 


2  2 

-Mr(x.y)  +  -Mf(x,y)  +  rf (xH(x,y)  =  0  , 

3x  3y^ 


(1.3.1a) 


where 


'i'(x.y)  =  I  A(x,y,p  ,p  T)exp{ii<()(x,y,p  ,p  )}dp  dp  , 
j  A  y  x  y  x  y 


(1.3.2a) 


or  alternatively, 


V2\p(r)  +  Tf(x)ip(r)  =  0  , 


(1.3.1b) 


where 


Mr)  =  A(r,p,T)exp(iT<ji(r,p)  }dp, 


(1.3.2b) 


—  —  —  —  — k 

A(x,y,px,p  ,t)  =  A(r,p,x)  -  E  Ak(r,p)t  =  T,  A^(x,y ,px,p^)t 

y  k  k 


i.e.,  r  =  (x,y)  ,  p=(pv.,p,). 

x  y 


■ 
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In  Chapter  II,  phase  space  analysis  is  summarized  and  the  Thom 
potentials  are  introduced.  In  Chapter  III,  <J> (r , p)  is  obtained  by  an 
extension  of  Maslov's  technique.  This  extension  allows  the  determination 
of  <}>(r,p)  for  media  characterized  by  a  wide  range  of  depth-dependent 
profiles.  Additionally,  this  approach  enables  explicit  representation 
of  the  emitter  geometry  and  the  boundary  conditions  in  cfi(r,p)  itself. 

The  geometry  of  the  caustic  in  configuration  space,  obtained  from  an 
analysis  of  the  stationary  points  of  4>(r,p)  in  phase  space,  is  also 
determined  in  Chapter  III.  In  Chapter  IV,  a  transport  equation  in 
mixed  configuration-momentum  space  is  derived.  This  transport  equation 
determines  the  higher  order  terms  in  the  asymptotic  series.  To  obtain 
the  asymptotic  series  of  the  integral  in  Equation  (1.3.2a),  alternatively 
(1.3.2b),  requires  a  coordinate  transformation  carrying  <J)(r,p)  to  an 
appropriate  canonical  form.  The  required  canonical  form  is  obtained 
from  a  criterion  based  on  Thom's  Classification  Theorem  (Appendix  A). 

In  Chapter  V,  the  coordinate  transformations  are  determined  explicitly. 

At  regular  points  in  the  field,  the  stationary  phase  technique 
determines  the  asymptotic  series  of  the  transformed  integrals.  At 
singular  stationary  points,  the  classical  stationary  phase  technique 
must  be  modified  to  determine  the  asymptotic  series.  In  Chapter  VI, 
the  asymptotic  series  at  both  regular  and  singular  stationary  points 
is  determined.  Chapter  VII  concludes  this  investigation  with  a  summary 
of  the  entire  algorithm,  an  explicit  example,  and  suggestions  for 


future  research. 
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HEURISTIC  SYNOPSIS  OF  THE  THOM  POTENTIALS 

2.1  Introduction 

Most  naturally  occurring  phenomena  are  nonlinear  in  character. 

Often,  an  exact  mathematical  analysis  of  such  phenomena  is  difficult; 
more  often,  it  is  only  possible  by  a  linear  approximation  of  the  non¬ 
linearity.  An  example  is  the  plane  pendulum,  where  the  nonlinear 
restoring  force  at  the  equilibrium  position  is  approximated  as  linear 
to  simplify  mathematical  analysis.  It  was  to  avoid  such  approximations 
that,  in  1881,  Poincare  [47]  introduced  the  qualitative,  or  topological, 
analysis  of  dynamical  systems. 

In  qualitative  dynamics,  a  dynamical  system  in  configuration  space 
is  represented  as  a  field  of  vectors  on  phase  space.  A  solution  is  a 
smooth  curve  tangent  at  each  point  to  the  vector  based  at  that  point. 
Qualitative  dynamics  is  concerned  with  the  topological  analysis  of  a 
family  of  solution  curves  throughout  the  entire  phase  space,  rather 
than  the  quantitative  analysis  of  a  particular  solution  curve  in  con¬ 
figuration  space.  The  family  of  solution  curves  in  phase  space  is  called 
the  phase  portrait. 

The  topological  analysis  of  a  phase  portrait  leads  to  a  global 
analysis  of  the  motion.  That  is,  the  entire  range  of  motion  of  the 
dynamical  system  can  be  studied  qualitatively  based  on  topological 
considerations.  Regions  corresponding  to  qualitatively  different 
motion,  e.g.,  periodic,  aperiodic  or  asymptotic,  are  easily  determined; 
and  the  separatrices,  or  divides,  which  demarcate  regions  of  stability 


and  instability,  are  readily  identified.  The  principal  limitation  of 
this  approach  is  that  it  does  not  lead  to  either  a  particular  explicit 
solution  curve  or  numerical  calculations.  Used  in  conjunction  with 
approximate  quantitative  methods,  such  as  asymptotic  expansions  in 
suitable  parameters  (also  developed  by  Poincare  [48]),  it  does  lead  to 
a  numerical  analysis  of  whatever  region  is  selected  for  study  [43]. 

The  phase  portrait  of  a  dynamical  system  is  obtained  from  the 
configuration  space  equation  of  motion.  The  equation  of  motion  is  re¬ 
duced  to  an  equivalent  set  of  first  order  differential  equations;  then 
an  integration  determines  a  solution  curve  in  phase  space.  The  family 
of  these  solution  curves,  each  member  corresponding  to  a  different 
initial  condition,  constitutes  the  phase  portrait.  For  example,  in  the 
second  order  differential  equation, 

x  =  F (x,x) 

where  the  dots  indicate  time  derivatives,  the  substitution  of  x  =  y, 
y=F(x,y)  leads  to 


An  integration  determines  the  solution  curve;  the  family  of  solution 
curves  is  the  phase  portrait,  which  is  represented  pictorially  as  the 
graph  of  x  vs.  x,  with  direction  indicated. 

2. 2  Structural  Stability  of  Differential  Equations 

In  1935,  Pontryagin  and  Andronov  extended  the  formulation  of 
Poincare  by  introducing  the  concept  of  structural  stability  [50].  A 
dynamical  system  is  said  to  be  structurally  stable  if , under  a  small 
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perturbation,  its  phase  portrait  is  qualitatively  unchanged.  For  example, 
if  a  dynamical  system  is  represented  by 


x  =  a^x  +  b^y  +  c^xy  +  . . . 
y  =  a^x  +  b^y  +  c2xy  +  . . . 


(2.2.1) 


the  system  is  structurally  stable  if  the  phase  portrait  is  qualitatively 
unchanged  under  a  small  variation  of  the  coefficients. 

Under  relatively  broad  assumptions  [13],  it  can  be  shown  that, 
in  a  neighborhood  of  the  origin,  the  analysis  of  Equations  (2.2.1)  can 
be  reduced  to  an  analysis  of 


= 

dx 


a2X  +  V 
a1x+b1y  ’ 


(2.2.2) 


Heuristically,  in  an  e-neighborhood  of  the  origin  (|ej  <  1), 
q11  >  qn+^  (n  >_  0) ,  where  q  is  either  x  or  y.  Hence,  the  linear  terms 
are  dominant  in  Equation  (2.2.1).  The  qualitative  nature  of  the  phase 
portrait  derives  from  the  eigenvalue  structure  of 


That  is,  the  eigenvalues  of  S  may  be  real,  imaginary  or  complex 
conjugates;  in  addition,  the  real  eigenvalues  may  have  the  same  or 
different  sign,  be  equal,  or  one  may  be  zero.  The  real  part  of  the 
complex  conjugates  may  be  positive  or  negative.  Each  eigenvalue 
structure  determines  a  qualitatively  different  phase  portrait  [43]. 

In  this  context,  strucutral  stability  refers  to  how  much  the  co¬ 
efficients  may  be  perturbed  without  altering  the  eigenvalue  structure 


of  S. 
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2.3  Structural  Stability  of  Potential  Functions 

Structural  stability  also  extends  to  phenomena  which  can  be 
characterized  by  a  potential.  In  this  context,  structural  stability 
is  closely  related  to  the  stability  of  equilibrium  points.  For  a  one¬ 
dimensional  system  described  by  a  potential,  V(q),  the  equilibrium 
point  is  determined  by  the  condition 

dvial  =  o 

dq 

d^V  , 

If  — —  f  0  at  the  equilibrium  point,  the  equilibrium  is  structurally 
dq 

stable;  however,  the  equilibrium  is  dynamically  stable  or  unstable  as 
d2V 

— —  is  less  than  or  greater  than  zero.  [In  this  investigation,  the 
dq 

historical  sign  convention  for  structural  stability  is  followed;  hence, 
the  force  is  defined  as  F =  V(q).}  If  at  the  equilibrium  position 

d2V 

— -r  =  0,  the  equilibrium  is  structurally  unstable.  That  is,  an 
dq 

arbitrarily  small  perturbation  of  the  potential  will  qualitatively  change 
the  phase  portrait  of  the  system. 

In  higher  dimensions,  the  equilibrium  point  is  determined  by 

VV(q)  =  0. 

The  equilibrium  is  stable  or  unstable  depending  on  the  eigenvalue 
structure  of  the  Hessian  matrix,  i.e.,  the  eigenvalues  of 
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at  the  equilibrium  point.  At  dynamically  stable  equilibrium  points, 
the  Hessian  determinant  is  zon-zero  and  the  eigenvalues  are  all  less 
than  zero.  If  the  Hessian  determinant  is  non-zero  and  the  eigenvalues 
are  not  all  less  than  zero,  the  equilibrium  is  dynamically  unstable; 
but  the  equilibrium  is  structurally  stable  in  that  an  arbitrarily  small 
perturbation  will  not  qualitatively  alter  the  phase  portrait  of  the 
system.  If  the  Hessian  determinant  is  zero,  i.e.,  the  Hessian  matrix 
has  at  least  one  vanishing  eigenvalue,  the  equilibrium  is  structurally 
unstable.  Such  a  point  is  called  degenerate  or  singular.  If  the 
Hessian  matrix  is  identically  zero  at  the  equilibrium  position  (all 
eigenvalues  zero),  the  equilibrium  is  completely  degenerate.  At 
structurally  stable  equilibrium  points,  it  is  possible  to  determine  a 
reversible  coordinate  transformation  carrying  V(q)  to  a  quadratic  form 
[44].  Such  a  coordinate  transformation  is  called  a  dif feomorphism. 

More  strictly,  a  dif feomorphism  is  a  differentiable  map  with  a 
differentiable  inverse.  Because  the  eigenvalue  structure  is  preserved 
under  dif feomorphism  [45],  the  quadratic  form  is  topologically  equivalent 
to  V(q).  The  condition  on  the  existence  of  the  dif feomorphism  is  that 


the  Hessian  determinant  is  non-zero  at  the  equilibrium  point,  i.e.,  that 
the  Jacobian  of  7V(q)  is  non-zero  at  the  equilibrium  point. 


2. 4  Heuristic  Synopsis  of  the  Thom  Potentials 

2.4.1  Background  and  introduction.  If  a  function  (potential)  with 
a  structurally  unstable  equilibrium  point  is  perturbed,  its  topology 
changes.  The  resulting  topology  is  not  always  determinable.  What  Thom 
did  was  to  determine  canonical  forms  for  a  class  of  functions  (those 
capable  of  parameterization  by  up  to  four  parameters,  e.g.,  space  and 
time)  which  could  be  perturbed  into  a  determinable  topology.  He  sought 
to  study  such  functions  by  including  them  as  isolated  members  of 
structurally  stable  families  and  analyzing  their  individual  behavior 
through  the  general  behavior  of  the  family.  Specifically,  Thom’s  Theorem 
(Appendix  A)  classifies  functions  into  one  of  seven  canonical  forms  based 
on  the  degree  of  degeneracy,  i.e.,  relative  structural  stability,  at 
the  equilibrium  point.  A  complete  treatment  of  the  theorem  requires 
rigorous  topological  analysis,  e.g.,  [18,53].  But  classical  dynamics 
does  provide  a  context  for  a  heuristic  synopsis,  as  noted  by  Haken  [29] . 


2.4.2  One-dimensional  potentials  (A-series).  Consider  the 
potential  V(q,A),  an  analytic  function  of  one  dependent  variable,  q, 
and  an  arbitrary  physical  parameter,  A.  For  a  given  Aq,  let  V(q,AQ) 
have  an  equilibrium  point  at  q =  0.  Then,  the  Taylor  series  of  V(q,AQ) 
about  the  equilibrium  point  becomes 


d2V(0) 


V(q,Ao)  =V(0,Aq)  +^-  + 


dq 


A(0)  3+...+i^^SW  %... 


jT - 


dq' 


n! 


dq 


The  equilibrium  is  dynamically  stable  or  unstable  as  the  sign  of 
2  2 

is  negative  or  positive.  If,  for  A  *  A  ,  — 
dq  dq^ 


=  0,  the  equilibrium 


17 


point  is  not  structurally  stable,  i.e.,  some  arbitrarily  small  pertur¬ 
bation  in  A  could  lead  to  either  a  stable  or  unstable  equilibrium. 

Let  the  parameter  A  be  such  that  at  A  =  A 

o 

dV (0)  =  d2V(0)  = 

dq  .2  U 

dq 

and 

^■#0. 

dq 

Then,  in  the  e-neighborhood  of  the  origin,  the  coordinate  transformation 

rl  d3V(0)  1  d4V(0)  1  dV0)  n-3  ,1/3 

5  q[3!  ,3  +4!  +-'-] 

dq  dq  dq 

transforms  V(q,AQ)  precisely  to 

V(5,Aq)  =V(0,Aq)  +C3  ,  (2. 4. 2.1) 

3 

with  the  coefficient  of  E,  normalized  to  unity.  (Heuristically ,  in 

an  c-neighborhood  of  the  origin,  q11  >>  qn+^;  hence,  the  potential  can 

be  closely  approximated  by  its  lowest  order  term.)  If  V(0,Aq)  =0, 

3 

V(C,Aq)=5  ;  this  is  the  simplest  single-variable  degeneracy. 

If  A  is  varied  slightly  to  A,  the  potential  becomes 


v(s,a)  =v(o,a) +ac  + e s2  +  £3  +  h4+  ...  , 


.3  . 


(2. 4. 2. 2.) 


where  the  coefficient  of  £  is  normalized  to  unity.  For  sufficiently 
small  A-variations  in  an  e-neighborhood  of  the  origin,  V(£,A)  may  be 
closely  approximated  by 
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V(5,X)  =  V(0,X)  +  a?  +ec2+YC3  ,  (2. 4. 2. 3) 

which  is  qualitatively  different  from  Equation  (2. 4. 2.1)  because  it  is 

nondegenerate,  i.e.,  structurally  stable,  at  its  equilibrium  point. 

Equation  (2. 4. 2. 3)  may  be  simplified  by  translating  the  origin  to 

2 

-8/3  and  shifting  the  zero-point  of  the  potential  to  28  /27  -aB/3  + 

V(0,X).  The  resulting  form  is 

V(£,X)  =  C3  +  <C  ,  (2. 4. 2. 4) 

2 

where  <  =  a  -  8  /9.  Equation  (2.4. 2.4)  is  the  form  of  the  lowest  order 

Thom  potential.  It  represents  a  structurally  stable  family  of 

potentials  parameterized  by  k.  The  structurally  unstable  potential  in 

Equation  (2.4.2. 1)  appears  as  an  isolated  member  when  <  =  0  and  the 

zero-point  is  appropriately  shifted.  [Thom  sought  to  study  the  behavior 

of  Equation  (2. 4. 2.1)  through  an  analysis  of  the  behavior  of  Equation 

(2.4. 2.4)  as  <  0. ]  Those  perturbations  which  qualitatively  change  the 

original  degeneracy  are  called  unfoldings.  Because  all  possible 
.  3 

unfoldings  of  £  near  the  origin  can  be  transformed  into  <£;, it  is 

called  the  universal  unfolding.  [While  this  treatment  has  been 

heuristic,  Thom's  Theorem  states  that  Equation  (2. 4. 2. 2)  may  be 

transformed  into  Equation  (2.4. 2.4)  exactly.] 

In  one  dimension,  there  exist  Thom  potentials  corresponding  to 

4  5  6 

q  ,  q  and  q  degeneracies  in  Equation  (2.4. 2.1).  An  analysis  similar 
to  the  above  determines  the  universal  unfoldings.  These  one-dimensional 
potentials  are  called  the  A-series  and  are  listed  in  Appendix  A. 
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2.4.3  Two-dimensional  potentials.  Consider  the  potential 
V(qi* q2»^) »  an  analytic  function  of  two  independent  variables,  and 
q^»  and  an  arbitrary  physical  parameter,  X.  For  a  given  X  =  Xq,  the 
equilibrium  point,  qQ  =  (q^Q  920^ »  is  determined  from 

VV(ql,q2’Ao)  =  0  • 

When  the  Hessian  determinant  is  non-zero  at  q^,  the  equilibrium  is 
structurally  stable.  If  the  Hessian  determinant  is  zero  at  qQ  (i.e., 
the  Hessian  has  at  least  one  vanishing  eigenvalue),  the  equilibrium 
is  not  structurally  stable.  When  the  Hessian  has  one  vanishing 
eigenvalue,  the  system  may  be  considered  to  lose  structural  stability 
in  one  direction  for,  under  a  principal  axis  transformation,  one 
homogeneous  quadratic  and  the  cross-term,  i.e,  9^2*  t*ie  Tayl°r 
expansion  of  V(q^,q2,XQ)  at  qQwill  vanish.  The  degree  of  degeneracy 
is  determined  by  the  degree  of  the  first  non-vanishing  Taylor  series 
term  in  the  direction  of  structural  instability. 

For  example,  let  V(q^,q2>XQ)  have  a  degeneracy  at  the  origin 
3 

corresponding  to  q  above.  Let  a(a  >  0,  for  clarity)  be  the  non-zero 
eigenvalue  of  the  Hessian  and  let  the  eigenvectors  be 


en" 

e2l" 

e  = 

e  * 

a 

_e!2_ 

o 

_e22_ 

Then,  the  principal  axis  transformation 
ql=  ?lell  +  P’2e12 


q2  =  ?le21  +  ^2e22 
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carries  VCq^q^A^  at  q  =  (0)  to 

V(C1,C2,Xo)  =  V(0,Ao)  +  g(C1,C2)C12  +  2;^22h(C2)  +  S-(C2)C23 
with  g(0)=a,  2.(0)  #  0.  Then,  the  coordinate  transformation 


e1  =  g(c1,c2)1/2c1  +  g(c1,c2)'1/2h(c2)c22 

-1  2  1/2 

B2=  U(C2) -g(?1,52)  h(C2)  S2]  /JC2  , 


where  the  negative  exponent  refers  to  powers  rather  than  inverses, 
determines  precisely 

V(61,B2,Ao)  =V(0,Aq)  +  612  +  623  . 

V(31>B2,Ao)  may  be  considered  the  canonical  form  for  the  degeneracy. 
The  universal  unfolding  in  the  direction  proceeds  as  in  the  one¬ 
dimensional  case.  The  analogous  canonical  form  is,  therefore 


V(B1,62,Aq)  =  B12  +  323  +  kB2  . 


When  the  Hessian  matrix  of  V(q^,q2,AQ)  at  q =  (0)  is  itself  zero, 
the  topology  of  the  phase  portrait  is  determined  by  the  cubic  terms  in 
the  Taylor  expansion  of  V(q^,q2>Ao)  about  q=  (0) 


V3(0,Ao)  =  c^3  +  c2qi2cl2+  c3qiq22  +  c4q23  ’ 


(2. 4. 3.1) 


where  V  (0,Aq)  represents  the  third  derivative  term  with  the  c^  the 
appropriate  Taylor  coefficients.  If  V  (0,Ao)#0,  and  at  least  one  of 
or  c^  is  non-zero,  the  cubic  is  equivalent  to  one  of  four  canonical 
forms.  The  equivalence  is  determined  by  an  analysis  of  the  root 
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3  — 

structure  of  the  cubic,  i.e.,  the  number  of  real  roots  of  V  (0,Xo)  =  0 
and  whether  or  not  any  of  the  roots  are  equal.  Each  canonical  form 
is  equivalent  to  a  Thom  potential. 

3  — 

To  determine  the  root  structure  of  Equation  (2. 4. 3.1),  V  (0,Xq) 

is  equated  to  zero,  then  reduced  to  a  one-dimensional  form  by  dividing 

3  3 

through  by  either  q^  or  q^  to  yield 

c^t3  +  c2t^  +  c^t  +  c^  =  0  ,  (2.4.3. 2) 

where  t  =  <1-^/ q2 ’  w*t*1  cl’c4^^  for  8enerality*  Equation  (2.4. 3. 2) 
can  have  four  root  structures:  three  real  equal  roots,  three  real 
roots  with  two  equal,  three  unequal  roots,  or  one  real  root  and  one 
complex  -  conjugate  pair.  Linear  transformations  have  been  determined 
which  carry  each  root  structure  to  the  canonical  forms  [49]: 


(i) 

three  real  equal  roots 

3 

:x 

(ii) 

three  real  unequal  roots 

3 

:x  -  xy 

(iii) 

one  real  root  and  one  complex 

3  3 

conjugate  pair 

:x  +  y 

(iv) 

three  real  roots,  with  two  equal 

2 

:x  y. 

x  is  the  lowest  order  Thom  potential.  When  associated  with  a  one- 

3 

dimensional  dynamical  system,  x  may  be  perturbed  into  the  unfolding 
3 

5  +  [Equation  (2. 4. 2. 2)].  When  associated  with  a  two-dimensional 

3 

dynamical  system,  a  perturbation  of  x  may  lead  to  any  of  an  infinite 

number  of  unfoldings,  i.e.,  an  indeterminate  topology.  Such  potentials 

are  qualitatively  different  from  those  considered  by  Thom  and  are 

3  2  3  3 

not  included  in  his  Classification  Theorem,  x  -  xy  and  x  +y 
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lead  to  the  unfoldings  which  are  called  the  elliptic  and  hyperbolic 

2  2  4 

unfoldings,  respectively,  x  y  leads  to  the  form  x  y  +  y  ,  which,  in 

turn,  leads  to  the  unfolding  called  the  parabolic  umbilic. 

If,  in  Equation  (2. 4. 3.1),  c^,c^ =  0  and  c^.c^ 5s 0,  the  cubic 

2 

corresponds  to  the  form  x  y.  If  c^.c^  an<^  one  c2  or  c3  ecIual  zero, 
there  is  no  equivalent  Thom  potential. 

3  2 

The  equilibrium  point  of  the  elliptic  umbilic,  cp  (x,y,A  )=x  -  xy  , 

L  O 

is  determined  from 

V(J>E(x,y,Ao)  =  0  , 

specifically 

(3x2  -  y2)i  -  2xyj  =  0  , 


to  be  (x,y)  =  (0).  At  (x,y)  =  (0),  the  Hessian  matrix  of  <J>  (x,y , A  )  is 

identically  zero;  hence,  it  is  completely  degenerate. 

If  A  is  varied  slightly,  e.g.,  A  -+A,  in  general,  *  (x,y,A  ) 

o  ho 

becomes  in  an  e-neighborhood  of  the  origin 


4>E(x,y,  A)  =  k1  + k2x  +  k3y  +  k^x2  +  k$y2  +  kfexy  +  x3  -  xy2  ,  (2. 4. 3. 3) 

3  2 

where  the  coefficient  of  x  -  xy  has  been  normalized  to  unity.  Then, 
by  translating  the  origin  to  x  =  x'  -k^/3,  y  =  y'  -k^/3  and  shifting 
the  zero-point  of  the  potential  to 


2  3  5 

kl  +  27  k4  +  9 


k4k5 

27 


k4k5 


+  k3k5, 


Equation  (2.4. 3. 3)  is  transformed  to 

3  2  2 

4>£(x,y, A)  =  x  -  xy  +K^y  +K2x  +  <3y  . 


(2. 4. 3. 4) 


Equation  (2. 4. 3. 4)  is  the  canonical  form  of  the  elliptic  umbilic. 

3  3 

An  analysis  of  the  hyperbolic  umbilic,  <{>_(x,y,X  )  =  x  +y  , 

proceeds  similarly.  The  equilibrium  point  is  again  (x,y)  =  (0);  at 

(x,y)  =  (0),  the  Hessian  matrix  of  <t>  (x,y,A  )  is  identically  2ero. 

H.  O 

If  X  is  varied  slightly  to  X,  in  an  e-neighborhood  of  the  origin, 
<t>H(x,y,Xo)  becomes 

2  2  3  3 

<J>H(x,y,X)  =  kx  +  k2x  +  k3y  +  k^x  +  k$y  +  k^xy  +  x  +y  , 


(2.4. 3. 5) 


3  3 

where  the  coefficient  of  x  +  y  has  been  normalized  to  unity.  Then, 
by  an  axis  translation  and  a  shift  of  the  zero  point  of  the  potential. 


Equation  (2. 4. 3. 5}  becomes 


3  3 

|>H(x,y,X)=x  +  y  +K1xy  +  K2x  +  <3y  , 


(2. 4. 3. 6) 


which  is  the  canonical  form  of  the  hyperbolic  umbilic. 

2 

The  form  (j>p(x,y,  X0)  =  x  y  cannot  be  treated  similarly,  for  applying 
the  equilibrium  condition  does  not  determine  a  unique  equilibrium  point. 
That  is,  setting 


yields 


V<J>p(x,y,Xo)  =  0 


2xyi  +  x  j  =  0 


At  x = 0,  the  condition  is  satisfied,  but  the  y-coordinate  is  not 
determined  uniquely.  Hence,  (J>p(x,y,X(j)  must  be  perturbed  into  an 
equivalent  form  which  determines  a  unique  equilibrium  point.  Because 
the  indeterminacy  is  in  the  y-coordinate,  a  perturbation  along  the 


24 


y-axis  seems  heuristically  appropriate.  Adding  +y  to  +y  to 
4>p(x,y,AQ)  uniquely  determines  the  equilibrium  point  but  also  obtains 

Hessian  matrices  which  are  not  completely  degenerate;  adding  +y  or 

2  3 

+y  renders  the  potential  structurally  stable.  Adding  +y  leads  to 

2  3 

x  y  +  y  ,  an  alternate  form  of  the  hyperbolic  umbilic  [18].  Adding 

3  2  3 

-y  leads  to  x  y  -  y  ,  the  form  of  the  elliptic  umbilic.  But  adding 
4 

y  to  4>p(x,y,Ao)  determines 

4>(x,y,Ao)  =x2y  +  y4  . 

At  (x,y)  =  (0),  the  Hessian  matrices  of  both  <(>p(x,y,Ao)  and  4>^(x,y,Ao) 

are  both  completely  degenerate;  however,  at  (x,y)  =  (0),  ^(x,y,Ao) 

has  a  unique  equilibrium.  In  an  e-neighborhood  of  the  origin, 

|x  y|  >  |y  |;  thus,  ^(x,y,AQ)  may  be  regarded  as  a  tractable 

2  4 

perturbation  of  4>^(x,y,Ao).  <Ji^(x,y, Aq)  =  x  y  +  y  leads  to  the  unfolding 

called  the  parabolic  umbilic.  (A  completely  rigorous  treatment  of  why 
2  ,24 

x  y  leads  to  the  form  x  y  +  y  requires  an  analysis  based  on  topological 
considerations,  e.g.  [18,35].)  Thus,  proceeding  as  with  the  elliptic 
and  hyperbolic  umbilics,  we  determine  the  form 


~  2  4  2  2 

t>p(x,y,A)  =  x  y +  y  +<1x  +*2y  +x3x  +  Ki!(y  , 


(2. 4. 3. 7) 


which  is  the  canonical  form  of  the  parabolic  umbilic. 

3 - 

If,  in  Equation  (2. 4. 3.1),  V  (q,AQ)  =0,  the  degeneracy  must  be 
studied  from  quartic,  or  higher,  terms  in  the  Taylor  expansion.  Such 
degeneracies  are  qualitatively  different  from  those  considered  above 


and  are  not  included  in  the  Thom  classification. 


This,  admittedly,  heuristic  synopsis  represents  an  amplification 
of  an  approach  outlined  by  Haken  [29].  The  results  of  this  chapter  have 
been  rigorously  proven  in  a  series  of  papers  by  Mather  [36-41]. 


CHAPTER  III 


DETERMINATION  OF  THE  PHASE 

3.1  Introduction 

The  determination  of  the  phase,  <j>(r,p),  and  its  relation  to  the 
medium  profile,  is  crucial  to  the  investigation  of  the  field  represented 
by  the  integral 

f 

4>(r)  =  A(r,p,T)exp{iTij>(r,p)  }dp  ,  (1.3.2b) 

4 

where 

A(r,p,t)  '  £  A  (r,p)f  k  . 
k  k 

In  the  mixed  configuration-momentum  space  approach,  the  phase  is  especially 
important.  Not  only  are  the  stationary  points  of  the  phase  obtained  from 
4,(r,p),  but  the  geometry  of  the  caustic  as  well.  The  standard  approach 
is  to  determine  <f>(r,p)  from  the  eikonal  Equation  (1.2.5)  and  the  A^(r,p) 
from  the  transport  Equation  (1.2.6)  in  momentum  space  [34].  As  in 
configuration  space,  in  momentum  space,  these  are  first  order  nonlinear 
partial  differential  equations. 

In  this  chapter,  the  phase  is  determined  for  a  medium  described 
by  a  single  variable  profile  through  an  extension  of  Maslov's  approach. 

This  procedure  will  allow  an  explicit  characterization  of  the  emitter 
geometry  and  the  boundary  conditions  in  the  phase  itself.  The  geometry 
of  the  caustic  in  configuration  space  is  then  determined  from  an 
analysis  of  the  stationary  points  of  <}>(r,p)  in  phase  space. 

3. 2  The  Technique  of  Maslov 

Maslov  considered  the  integral  representation  for  the  field, 

Equation  (1.3.2b).  But  rather  than  determine  the  phase  from  the 


eikonal,  he  showed  <(>(r,p)  had  the  form 


<t>(r,p)  =  r  •  p  -  S (p)  , 
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(3.2, 


where  S(p)  is  the  generating  function  of  a  canonical  transformation  [3] 
From  Equation  (3.2.1),  the  stationary  phase  condition,  V^4) (r , p)  =  0, 
becomes 


r-  VpS(p)  =  0, 


(3.2, 


i.e.,  a  Lagrange  manifold.  A  Lagrange  manifold  may  be  defined  as  a 
surface  determined  by  the  gradient  of  a  generating  function.  By  re¬ 
writing  Equation  (3.2.2)  in  component  form,  i.e., 


_  as(P)  = 


v  _  9,s(P.) 

y  3o 

•  y 


o 


(3.2 


*  0  , 


the  Lagrange  manifold  may  be  seen  as  a  transformation  from  momentum 
space  to  configuration  space.  On  the  Lagrange  manifold,  i.e.,  for 
r  =  V^S(p),  the  eikonal  condition  in  Equation  (1.2.5)  becomes 


p  •  p  -  f  (r)  *  0  , 


or  more  explicitly 


p  •  p-f[VpS(p)]  =0  . 


Maslov  proceeded  by  defining  his  Hamiltonian  as 


(3.2 


(3.2 


H  =  p  •  P  -  f (r)  , 


(3.2 


i) 


2) 


•  3) 


.4a) 


.4b) 


•5) 


for  arbitrary  (r,p).  Because  the  Lagrange  manifold  is  invariant 
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under  the  Hamiltonian  flow  determined  from  Equation  (3.2.5),  i.e.,  a 
Hamiltonian  trajectory  satisfying  Equation  (3.2.2)  for  t  =  0  satisfies 
Equation  (3.2.2)  for  all  t  [10],  S(p)  may  always  be  determined,  e.g.,  using 
Hamilton-Jacobi  theory  [31]. 

For  example,  in  two  dimensions  S(p)  may  be  determined  by  applying 
Hamilton's  equations  to  Equation  (3.2.5),  leading  to 

=  5f (*>y) 


x  =  2p 


y  =  2p 


'  x  8x 

’  =  (x,y) 

y  3y 


which  obtains  the  map 


x  “  x(t ,  0)  p  =p  (t,0) 

X  X 


y  =  y(t,0)  py  =  py(t,e)  , 


(3.2.6) 


where  t  is  the  time  and  0  is  a  parametrized  initial  condition.  Then 
inverting  the  momentum  equations  in  Equation  (3.2.6)  yields 

(3.2.7) 

Q  =  o(px>py). 


Substituting  Equation  (3.2.7)  into  the  configuration  space  equations 
in  Equation  (3.2.6)  determines 

x  =  x[t(p  ,p  ),  0(p  ,p  )]  =X(p  ,p  ) 

X  y  X  7  X  y  (3.2.8) 

y  =  y[t(px,py),  0(px,py)l  =  Y(px,py)  , 

where  X(px»py)  and  Y(Px*Py)  are  explicit  functions  of  p^  and  p^. 

Since  the  Hamiltonian  flow  preserves  the  Lagrange  manifold,  from 
Equations  (3.2.3)  and  (3.2.8), 
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m 


iu  ip  1 1 


^  -  *<W 


■  Y(pk’V  • 

y 


(3.2.9) 


Integrating  Equations  (3.2.9)  determines  S(p)  which,  in  turn,  determines 
4>(r,p),  [28]. 

3.3  An  Extension  of  Maslov's  Technique 

For  those  media  characterized  by  a  single  variable  profile,  f(x), 
an  extension  of  this  approach  leads  to  a  phase  which  includes  an 
explicit  representation  of  the  emitter  geometry  and  the  boundary 
conditions  on  the  wave  vectors. 

Consider  the  Hamiltonian  (on  the  Lagrange  manifold) 


H=p  2  +  p  ^-f (x)  =  0  . 
x  *y 


(3.3.1) 


Then,  at  any  field  point  where  f  ^  is  defined, 


-12  2 
x  =  f  (Px  +  Py  ) 


(3.3.2) 


determines  one  coordinate  of  the  Lagrange  manifold 

x=  igill 

3px  ' 


(3.3.3) 


The  other  spatial  coordinate  may  be  determined  by  noting  that,  for 
2 

gradient  (C  )  functions, 

52S(p)  _  32S  (p)  . 

3VPy  3py3px  ’ 


hence, 


d  _  +  0(p  )  =  3S(J 2l 

3Px3Py  Px  3Py  °(Py)  3Py  ’ 


(3.3.4) 
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where  G(p^)  is  an  arbitrary  function  of  p^.  Thus, 


3p 


SfiEl  ,  +  e 

3py  »py  y 

define  the  Lagrange  manifold,  where  G(p^)  derives  from  the  boundary 
conditions  at  the  emitter. 

Let  ©(Py)  be  represented  by  a  polynomial  in  p^  of  arbitrary 
power  with  constant  coefficients,  and  let  the  geometry  of  the  emitter 
be  represented  by 

y  =  g(x)  .  (3.3.5) 


The  spatial  coordinates  of  the  emitter  taken  together  with  the  initial 
boundary  conditions  on  the  wave  vectors  may  be  regarded  as  a  lifting 
from  configuration  space  to  the  mixed  configuration-momentum  space, 
i.e.,  from  2-space  to  4-space.  Then,  by  choosing  the  intital  conditions 
on  the  emitter  so  that,  at  t =  0,  [r(t),p(t>]  satisfy  Equation  (3.3.2), 
Equations  (3.3.3),  (3.3.4)  and  (3.3.5)  determine 

g(x),lL(fI+  a  +a  p  +  .  as (I). 

8W  3p  +ao  +  alPy+ 3p  ’ 


(3.3.6) 


on  the  emitter.  Solving  Equation  (3.3.3)  for  px  and  substituting  into 
Equation  (3.3.6)  determines  an  equation  in  x  and  p^.  Successive 
differentiations  with  respect  to  x  lead  to  a  system  of  linear  algebraic 
equations  for  the  a^  in  terms  of  the  initial  wave  vectors  and  the 
spatial  variations  of  the  wave  vectors  at  a  point  on  the  emitter. 
Solving  for  the  a^  determines  Equation  (3.3.4)  explicitly.  Equations 
(3.3.3)  and  (3.3.4)  determine  the  Lagrange  manifold,  and  hence, 


T 


IWM^  ^  ■WjRp-k' 
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’(r,p)  =  4>(x,y,px,py)  =  xpx  +  ypy  -  S(px,py)  . 


(3.3.7) 


3.4  Determination  of  the  Caustic  Geometry 


The  analysis  of  the  field  at  any  point  rQ =  (xQ,yQ)  also  proceeds 

from  the  phase.  The  Lagrange  manifold  conditions  at  rQ  are  solved 

for  the  associated  momentum,  p  =  (p  ,p  ).  If  the  Hessian  determinant 

o  xo  yo 

of  cf>(r,p)  at  (ro,pQ),  i.e., 


det 


a2d> 

32d> 

3p2 

*x 

apx3Py 

9 

2 

3  d> 

3d> 

3VPx 

3  2 

py 

-<VPo) 


is  non-zero,  (r^p^)  is  a  regular  point.  The  non-vanishing  of  the 
Hessian  determinant  at  a  stationary  point  is  analogous  to  the  non¬ 
vanishing  of  the  second  derivative  of  the  phase  in  one  dimension 
(Chapter  I).  Hence,  it  is  the  condition  for  the  validity  of  the  classical 
stationary  phase  technique  in  higher  dimensions  [15].  If  the  Hessian 
determinant  of  <>(r,p)  at  (r0»PQ)  is  zero,  i.e.,  <p( r,p)  is  structurally 
unstable  at  (r0’PQ)>  (r^.p^)  a  s:*-nSu^ar  P°int*  Consequently,  the 
classical  stationary  phase  technique  would  not  be  valid  at  (r^p^). 

The  equation  of  the  caustic  may  then  be  determined  from  the 
Hessian.  Equating  the  Hessian  determinant  to  zero  defines  a  curve  of 
singular  points  in  momentum  space,  i.e.,  the  caustic  in  momentum 
space.  Associated  with  every  point  p^  on  the  caustic  is  a  point 
ro*(xo,yQ)  in  configuration  space  determined  by  the  Lagrange 


manifold  condition 
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x 

o 


3S(Pq) 

3px 


3S(pq) 

"3PT” 


The  locus  of  these  points  determines  the  equation  of  the  caustic  in 
configuration  space. 


v  CHAPTER  IV 
THE  TRANSPORT  EQUATION 

4.1  Introduction 

In  this  investigation*  the  reduced  Helmholtz  wave  equation 

v2<Kr)  +  Tf(x><Hr)  =  Q  ,  "  (1.3. lb) 

where  t  is  a  large  parameter,  has  an  asymptotic  solution  of  the 
form 

yi(r)  —  A(r ,p,x)exp{ it<j>(r ,p)  }dp  =  0(x_°°)  ,  (4.1.1) 

where 

A(r,p,x)  «  l  A  (r,p)x  k 
k  k 

which  is  an  alternate  form  of  Equation  (1.3.2b),  [22,23].  Ordinarily, 
the  phase  is  determined  from  the  eikonal  Equation  (1.2.5)  and  the 
series  coefficients  are  determined  from  the  transport  Equation  (1.2.6). 
Because  the  eikonal  and  transport  equations  have  representations  in 
both  configuration  space  and  momentum  space,  it  is  possible  to  determine 
the  phase  and  the  series  coefficients  exclusively  in  either  space  [24]. 
Here,  following  Maslov  and  Duistermaat  [22,23],  we  consider  a  mixed 
configuration-momentum  space  representation  for  the  phase,  <j>(r,p),  and 
the  amplitude,  A(r,p,x).  The  phase  has  been  determined  without  recourse 
to  the  eikonal  equation  (Chapter  III).  To  determine  the  series 
coefficients,  A^(r,p),  it  is  necessary  to  derive  an  equation  in  the 
mixed  space  representation  that  is  analogous  to  the  transport 
equation. 
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4.2  Derivation  of  the  Modified  Transport  Equation 

Let  the  asymptotic  solution  of  Equation  (4.1.1)  be  represented 


as 


*(r)  - 


A(r ,p,T>exp{ It [r  •  p  -  S(p)]}dp  =  0(t  ) 


(4.2.1) 


i.e. ,  the  form  of  the  phase  (from  Chapter  III)  is  explicitly  noted. 
The  derivation  of  the  transport  equation  proceeds  by  carrying  the 
dif f erentation  in  Equation  (1.3.1b)  across  the  integral  in  Equation 
(4.2.1),  obtaining 

V2iJ;(r)  +xf  (x)!|;(r) 


exp{  ix  [r  •  p-  S  ( p)  ]  }  f  (  ±t  ) 2  (p  •  p-  f(x))A+2ix(p  •  V  A) 

r 


+  Vr2A]dp  =  0<t"")  , 


(4.2.2) 


where,  for  clarity,  the  argument  of  A(r,p,x)  is  included  implicitly. 
Expanding  p  •  p  -  f (x)  in  a  Taylor  series  about  the  turning  point 
determines 


p  •  p  -  f  (x)  =  p  •  p  -  f(3j>xS)  +  {r- V  S(p)]  •  D(r,p)  ,  (4.2.3) 


where  D(r,p)  is  the  remainder  term,  which  may  be  represented  as 

rl 


D(r,p)  =  D  =>  - 


V  f [t(r-V  S) + V  S ]dt  . 
r  P  P 


The  first  term  is  precisely  Maslov's  Hamiltonian  [and  thus  is  zero 
by  Equation  (3.2.6)]  on  the  Lagrange  manifold;  hence.  Equation  (4.2.3) 


becomes 
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p  •  p  -  £ (x)  =  [r  -  VpS(p)]  •  D(r,p)  .  (4.2.4) 

Substituting  Equation  (4.2.4)  into  Equation  (4.2.2)  determines 

V2^(r)  +  tf  (x)4»(r) 

* 

-  exp{iT[r  •  p  -  S(p) ] } [ (it)2 [r  -  VpS]  •  DA+  2it(p  •  VrA) 

+  Vr2A]dF  =  0(t-°°)  .  (4.2.5) 

But 

exp{ it [r  •  p  -  S(p) ] } (it) 2 [r  -  7pS(p) ]  •  DAdp 
=  it  •  (exp{it[r  •  p  -  S(p) ] }AD)dp 

4 

' 

-it  exp{ir[r  •  p  -  S(p)]}[D  •  V^A  +  AVp  •  D]dp  .  (4.2.6) 

By  substituting  Equation  (4.2.6)  into  Equation  (4.2.5)  and  taking  the 
surface  integral  over  a  sufficiently  large  radius  that  it  vanishes, 
Equation  (4.2.5)  becomes 

V2\Jj(r)  +  rf  (x)tp(r) 

-  exp{ it [r  •  p  -  S(p) ] }it [-D  •  V  A  -  AV  •  D  +  2p  •  V  A 

P  P  r 

. 

+  ^  Vf2A]dF  =  0(t-°°)  .  (4.2.7) 

For  Equation  (4.2.1)  to  have  an  asymptotic  solution,  it  is 
necessary  that  Equation  (4.2.7)  be  satisfied  in  a  neighborhood  of 
the  Lagrange  manifold.  Noting  Equation  (1.3.1b),  a  sufficient  condition 
for  Equation  (4.2.7)  is  that 


36 


-D  • 


V  A-AV  •  D+  2p  •  V  k+~  V  2A  =  0. 
p  p  r  it  r 


(4.2.8) 


Then  by  introducing  the  flow 
r  =  2? 

p  =  -D(r,p) , 

Equation  (4.2.8)  becomes 

V\V5+'A-l-°  ' 


(4.2.9) 


(4.2.10) 


where  the  dots  indicate  the  time  derivatives  determined  by  the  flow 
in  Equation  (4.2.9).  Equation  (4.2.10)  determines  the  evolution  of 
the  A^'s  as  does  the  transport  equation  (in  either  configuration 
or  momentum  space  [24]).  Hence,  we  regard  Equation  (4.2.10)  as  the 
transport  equation  in  the  mixed  configuration-momentum  space 
representation.  [It  should  be  noted  that  if  the  flow  is  defined  by 
the  Hamiltonian,  i.e., 

1=2? 

P  =  Vxf(x)  , 

it  is  only  possible  to  determine  the  asymptotic  limit,  i.e.,  A  (r,p).] 

4.3  Outline  of  the  Algorithm 

The  basic  algorithm  can  now  be  outlined.  At  any  point,  rQ =  (xQ,yo), 
in  the  field  the  phase  is  determined  by  an  extension  of  Maslov's 
technique  (Chapter  III).  Next,  AQ(r,p)  is  obtained  from  the  transport 
equation.  Then,  the  asymptotic  series  of  the  integral 

Ao(r0,P)  exp{it<}i(ro,p)  Jdp 
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is  determined.  Additional  A^(rQ,p)  proceed  from  Equation  (4.2.10), 
yielding  the  integrals 

Ak(rQ,p)exp{iTij)(ro,p)}dp  . 

. 

Each  integral  determines  an  asymptotic  series.  The  entire  asymptotic 
series  at  rQ  is  the  sum  of  the  series  of  the  integrals. 


I 


CHAPTER  V 


THE  COORDINATE  TRANSFORMATIONS 

5.1  Introduction 

We  consider  an  integral  of  the  form 

1=  A(x,y,p  ,p  ,T)exp{iT<j>(x,y,p  ,p  )}dp  dp  ,  (1.3.2a) 

%  x  y  x  y  x  y 

. 

where  A(x,y,p  ,p  ,t)  and  all  its  derivatives  are  bounded  throughout 
x  y 

the  region  of  integration,  A(x,y,px,p^,T)  is  analytic  near  the  stationary 

points  of  (j)(x,y,p  ,p  )  and  (j)(x,y,p  ,p  )  has  no  poles  in  the  region  of 
x  y  x  y 

integration.  Before  the  asymptotic  series  of  Equation  (1.3.2a)  can  be 

determined,  the  integral  must  be  transformed  to  a  form  appropriate 

for  asymptotic  analysis,  i.e.,  the  phase  must  be  transformed  to  a 

canonical  form.  The  transformation  proceeds  on  the  nature  of  the 

Hessian  determinant  of  the  phase,  <j>(x,y,p  ,P  ),  at  the  stationary 

x  y 

point,  (r  ,p  )  =  (x  ,y  ,p  ,p  ).  If  the  Hessian  determinant  is  non- 
r  o  ro  o^o  xo  yo 

zero  at  the  stationary  point,  <Ji(x,y,p  ,p  )  at  (r  ,p  )  is  transformed 

x  y  o  o 

to 

^ (Fo,  B1,  S2)  =  <j>(Fo,F0)  ±  Bx2  ±  622  ,  (5.1.1) 

which  follows  from  the  Morse  Lemma  [44].  If  the  Hessian  determinant 
is  zero  at  the  stationary  point  with  one  vanishing  eigenvalue,  then 
the  phase  is  transformed  to 

=<!>(^0,P0)±8l2±32n  »  (5.1.2) 

which  follows  from  the  Gromoll-Meyer  Lemma  [27].  The  exponent  of  0 2 
is  determined  from  the  Thom  potential  to  which  <}>(r,p)  at  (rQ,Po) 
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is  equivalent  (Appendix  A).  When  the  Hessian  determinant  at  the  station¬ 
ary  point  has  two  vanishing  eigenvalues,  the  transformed  potential  assumes 
a  qualitatively  different  form  than  in  Equations  (5.1.1)  and  (5.1.2). 

This  case  is  not  considered  here.  The  transformation  of  the  phase  to  a 
canonical  form  allows  the  integral  in  Equation  (1.3.2a)  to  be  expressed  in 
a  form  suitable  for  asymptotic  analysis.  In  this  chapter,  we  determine 
the  required  coordinate  transformations. 


5.2  Regular  Stationary  Points 

At  rQ  =  (xo,yQ),  let  <j>  (x.y,Px»Py)  =  4>  (r , p)  have  a  stationary  point 

[p  =  (p  ,p  )  such  that  V  $(r  ,p  )  = 0,  taken  here  at  p  =  (0)  for 
o  xo  yo  p  o  o  *o 

clarity].  Then,  the  Taylor  series  of  q> (r , p )  at  (r0»PQ)  is 

4>(r,p)  =  <j)(ro,0) +C1(px)px2  +  2pxC2(px,py) +C3(py)py2  ,  (5.2.1) 


where  <^(0)  and  0^(0)  ^  0.  Let 


®Xl 


"11 


'12 


*X2  = 


'21 


'22 


be  the  eigenvectors  corresponding  to  the  eigenvalues  at  X^  and  X^, 
respectively.  Then,  the  principal  axis  transformation 


(5.2.2) 

Py  =  ?le21  +  S2e22 

and  Taylor's  Theorem  carry  (r , p)  at  (r0*P0)  to 

“  4>(ro,0)  +g1(C1,C2)C12  +  2C1f,22h(?2)  +g2(C2)C22  , 
where  g^O)  and  g2(0)  are  the  eigenvalues.  Then,  by  completing  the 


square,  we  determine  the  coordinate  transformation 


AO 


1/2  -1/2 

S1  =  ^1(C1,C2)  ?i±8i^i»52)  hU2>V 

-1  2  1/2 

e2 =  [g2<C2)  +8i(?i.C2)  h(c2)  ]  C2  » 


(5.2.3) 


where  g.  (?)  =  g-j  (O  when  g-L(0)  >0,  and  g1(?)=-g1(0  when  g, (0)<  0.  The 
positive  sign  in  g^  occurs  when  the  Hessian  of  <}>(r,p)  is  positive 
definite  at  (r  , p  ).  The  positive  sign  in  g„  occurs  when  the  Hessian 
of  4> (r , p )  is  definite  indefinite  at  (r0»P0)»  with  g2(0)  the  positive 
eigenvalue.  Otherwise,  the  negative  signs  appear  in  the  transformation. 

At  (r0»P0)»  Equation  (5.2.3)  carries  <f>(r,p)  at  (r0’PQ)  t0 

^(7q,61,B2)  =  <|>(ro,0)  B22  ,  (5.2.4) 

where  the  sign  of  each  g^  is  determined  by  whether  the  corresponding 
eigenvalue  is  positive  or  negative. 


5.3  Singular  Stationary  Points 

If  the  Hessian  determinant  of  <j>(r,p)  is  zero  at  the  stationary 

point,  <(>(x,y,p  »P  )  at  (r  ,p  )  is  transformed  to 
x  y  o  o 

$(ro,g1,g2)  =  <t>(ro,o)  ±  e12±  g2n  , 

where  n>^3.  If  n  =  3  or  4,  the  coordinate  transformation  carrying 
4>(r,p)  at  (r0>PQ)  to  Mr0>3^»62)  may  again  be  obtained  by  explicit 
algebraic  computation.  For  n>4,  the  determination  of  the  appropriate 
coordinate  transformation  requires  the  algorithm  specified  in  the 
proof  of  the  Gromoll-Meyer  Lemma,  i.e.,  implicit  rather  than  explicit 
solutions  of  equations,  e.g.,  [49].  The  procedure  can  be  made 
explicit,  however,  by  an  application  of  the  Cauchy  inversion 
formula  [36]. 
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Let 


"iT 

1 - 

rH 

CN 

0) 

L  _ 

_e!2  _ 

* 

n> 

o 

u 

_e22  _ 

be  the  eigenvectors  corresponding  to  the  eigenvalues  at  X  and  0  (zero) , 
respectively.  Then,  the  principal  axis  transformation  in  Equation  (5.2.2) 
and  Taylor's  theorem  carry  <j>(r,p)  at  (rQ,po)  to 


*(  VW  ~  *<v0)  +  ^i»s2Ki2  +  2^22mv  +  ns2n2n,  (5.3.D 

where  n  equals  3  or  4,  g(0)  =X  and  it (0)  ^  0 .  When  n=3,  by  completing 

the  square,  we  determine  the  coordinate  transformation 

1/2  -1/2 

B1“g(C1.52)  51±i(C1,C2)  h(?2)52 

-1  1/3  (5.3.2) 

ha2)  *2]  K2  ’ 


where  g(C1, C2)  =  g(C1» ?2)  when  g(0)  >  0,  and  gU^ ?2)  =  -g(C1, when 

g(0)  <0.  The  sign  in  0^  is  determined  by  the  sign  of  X.  At  (rQ,po), 

Equations  (5.3.2)  carry  4> (r  ,p  )  to 

o  o 

MFo,B1,e2)  =  <KFo,0)  +  B12+623  ,  (5.3.3) 

where  the  sign  of  0^  is  determined  by  sign  of  X.  Similarly,  when  n=4, 

we  determine  the  coordinate  transformations 

1/2  _  -1/2  , 

e1  =  g(c1,c2)  h(c2)e2 

(5.3.4) 

-1  2 

62=  [£(52) -g(51,Cz)  h(C2)  ]  C2  , 

for  Jt(0)  -  g(0>_1h(0)2  >  0  and 
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1/2  -1/2  , 
e1  =  g(c1,€2)  c1±i(c1,e2)  h(c2)c2 


(5.3.5) 

-1  2 

B2— IgUj.Cj)  h(C2)  -2(C2)]  ?2 

for  S,(0)-g(0)  1h(0)2<0,  where  g(C1>C2)  =  g(C1>?2>  when  g(0)  >  0,  and 
g(^l’^2^  *  “S^i’^)  w^en  8(0)  <  0*  The  positive  sign  in  6  corresponds 
to  X  >  0.  At  (r0*P0)»  Equations  (5.3.4)  or  (5.3.5)  carry  <j>(ro,Po)  to 

t(ro,61,02)  =<j>(7o,O)  +  012  +  g2A  ,  (5.3.6) 

where  the  sign  of  g^  is  determined  by  sign  of  X  and  the  sign  of  g 

- 1  2 

is  determined  by  the  sign  of  £(0)  -  g(0)  h(0)  . 

When  n>4,  explicit  algebraic  computation  does  not,  in  general, 
suffice  to  determine  a  coordinate  transformation  which  carries  <j>(r,p) 

at  (VPo}  t0 


61,g2)  =  <Kro,0)  +  Bj2  +  3kn  . 


The  procedure  for  determining  the  transformation  can  be  made  explicit, 
however.  Application  of  the  principal  axis  transformation  in  Equation 
(5.2.2)  to  <j>(r,p)  at  (r0»PQ)  determines 

M7q,C1,C2)  =  4>(ro,0)  +  g(C1,C2)C12  +  2S1C22h(C1,C2)  +  JK£2K2n,  (5.3.7) 


where  g(0)  is  the  non-vanishing  eigenvalue  and  £(0)  0.  Because 

2  —  —  2 - 

3^(0)  1*0,  i.e.,  g(0)  i*  0,  3^  <P(0)  i*  0;  hence,  from  the  Implicit  Function 
Theorem  [21],  setting  3^4i(C^»C2)  =  0  determines  an  implicit  equation 
for  in  terms  of  £2»  i.e., 

c1  =  o(?2)  • 
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Q (C 2)  may  be  determined  explicitly,  however,  from  the  Cauchy  inversion 
formula  [35],  i.e.. 


e(«2>-2H 


z312ij!(z,C2)dz 
3^(2, C2) 


(5.3.8) 


for  fixed  E,^,  where  y  is  a  circle  in  C2  =  0  enclosing  =  0.  Then 
introducing  the  transformation 


5i  =  ai  +  0(a2) 


?2  a2  > 


Equation  (5.3.7)  becomes 


(5.3.9) 


<t>[ro’ai  +  0(a2^’ail  =  4>  (r0» +  0(a2)  ,a2l  [o^+G^)]  + 

2[a1  +  0(a2)]a22h[a1  +  0( a2),«2]  +  £(a2)a2n  . 


(5.3.10) 


Because  the  Taylor  series  of  an  analytic  function  of  two  variables  may 
be  written  as  a  Taylor  series  in  one  variable  with  the  other  fixed  [35], 
expanding  Equation  (5.3.10)  about  =  0  with  a2  fixed  determines 


4>[ro.«i  +  0(«2)  >a2]  =  4>(ro,0)  +  $[ro,0(a2)  ,a2]  + 

-  2 

al3l4>[ro,0^a2^  ,a2^  +al  R(ai»a2^  » 


(5.3.11) 


where  R(a^,a  )  is  the  remainder  of  the  Taylor  series  less  a  factor  of 
2  —  2  —  — 

,  and  R(0)  3*  0  [since  3^  <J>(0)/0],  Because  9^<{>(C^,C2)  =  0,  Equation 
(5.3.11)  becomes 

4>[r0»ai  +  0(<*2)  »«2^  =  (r0 » 0)  +  4>[rQ,0(u2)  ,«2]  +  a12R(a1,a2)  .  (5.3.12) 

Then  the  transformation 
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Yx = a1|R(a1,a2) 


1/2 


(5.3.13) 


Y2  =  Ct2 


obtains 


<Kro,Y^»Y2)  =  <Kro»0)  i^i  +  u(y2>  =  +  0(a2)  ,a2]  . 


(5.3.14) 


Finally,  since  the  first  non-vanishing  Taylor  coefficient  (at  ?2 = 0) 
in  the  ^-direction  is  the  n-th  [Equation  (5.3.7)],  u(y2)  may  be 
expressed  as 


u(Y2)  = Y2ny(Y2)  » 


where  p(0)^0.  Consequently,  the  transformation 


B1  =  Y1 


B2  =  Y2'm(y2^ 


1/n 


(5.3.15) 


carries  <Kr  ,y.  ,yJ  to 
o  i  z 


^(r0.Bi.32)  =  <Kro.0)±B12  +  82n  » 


(5.3.16) 


where  the  signs  of  the  0^  must  be  determined  for  each  specific  case, 
as  above.  Equation  (5.3.16)  is  the  form  of  the  phase  required  for 
the  determination  of  the  asymptotic  series,  c.f.  Equation  (5.1.2). 

5.4  Transformation  of  the  Amplitude 

Under  the  coordinate  transformation  which  carries  $(r,p)  to 
the  appropriate  canonical  form,  the  integral  in  Equation  (5.1.1)  at 
(rQ,po)  becomes 


i  '  ; 
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A(x,y,p  , p  , x)exp{ (x,y,p  ,p  ) }dp  dp 
x  y  x  y  x  y 


(5.4.1) 


2  .  „ 


exp{iT<}i(ro,po)}  A(ro,$1,32,x)exp{iT(+S1  +  B2  )}dB.jd82  , 


_  3(p  ,p  ) 

where  A(rQ,  8^  B2»t)  =  A(rQ,px>py,T)  •  »  which  is  the  form 

required  for  the  determination  of  the  asymptotic  series.  Since  A(r,p,t) 
is  analytic  over  the  region  of  integration,  ACr^B^^T)  may  be 


expressed  as 


i(ro>6i,62,t)  -  Z  brs6l  S2  _ 


i  ffA(v>wT)'ifylpJ 


where 


rs  .2  t  Nr+1  /  vS+1  * 

^  JJVvV  y2(px’V 


Yl“0l[5l(p*’V’  C2(px’py)] 


Y2=¥5l(px’V’  C2(px’py)]  , 


(5.4.2) 


(5.4.3) 


(5.4.4) 


and  6  and  e  are  contours  enclosing  the  origin.  Equation  (5.4.3),  which 

determines  the  b  explicitly,  proceeds  from  an  extension  of  the  Cauchy 
r  s 

inversion  formula  to  two  variables. 

More  specifically,  in  one  complex  variable,  if  f(z)  is  analytic 
and  single-valued  at  zq,  and  f  ^(z)  exists  and  is  single-valued  in  a 
neighborhood  of  w  =f(zQ),  the  Cauchy  inversion  formula  takes  the 


form 
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g(f  1(w)) 


1  g(z)f 1 (z)dz 
2iri  f(z)-w 

° 

Y 


(5.4.5) 


where  y  is  a  circle  enclosing  zq.  If  is  an  analytic  function 

of  two  complex  variables,  it  may  be  regarded  as  an  analytic  function 
of  z^  for  fixed  z 2  as  an  analytic  function  of  z 2  for  fixed  z^  [35]. 

The  Implicit  Function  Theorem  implies  that  two  analytic  functions  each 
of  two  complex  variables,  e.g.. 


W1  =  81^Z1’Z2^ 


w2  =  g2(z1,z2) 


with 


3(w1,w2) 

3(z1»z2) 


(0)  *0 


9 


may  be  inverted  to  determine 
z1  =  ©1(w1,w2) 

(5.4.7) 

z2  =  02(w1,w2)  . 


By  repeated  application  of  the  Cauchy  inversion  formula,  however  Equation 
(5.4.7)  may  be  determined  explicitly.  Heuristically,  first  z^  is 
determined  as  a  function  of  z2  and  w^,  e.g.,  z^*h(w^,z2).  Next  z2 
is  determined  as  a  function  of  w^  and  w^.  Finally,  substituting  for 
z2  in  z^  =  h(w^,z2)  determines  z^  as  a  function  of  w^  and  w2* 

More  explicitly. 


Z1  =  h^wi‘z^ 


1 

2vi 


n31g1(n,z)dn 

g1(n,z)-w1 


c 


(5.4.8) 
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where  c  is  a  contour  enclosing  the  origin  and,  for  clarity,  it  is 
assumed  g^O)  =  g2(0 )  =  0  and  3^(0)  i  0.  Since  w 2  =  g2^zl,z2^ 

w2 = g2[h(w1,z2),z2]  ,  (5.4.9) 


z0  =eo(w .  ,w„)  = 


f{31g2[h(w1,n) »n]a2h(w1,n)  +  a2g2[h(w1>n) ,n]}ndn 


2  2  1  2  2-ni 


g2lh(w  ,n) , n]  -  w2 


(5.4.10) 


where  c'  is  a  contour  enclosing  the  origin.  Then 


z1 = 0(w1,w2)  =  h[w1,02(w1,w  ) ]  . 


(5.4.11) 


An  analytic  function  of  two  complex  variables,  f(z^,z2),  may  be 
written  directly  as  a  function  of  two  other  complex  variables, f(w^,w2) , 
x.  e. , 

j 

f  (z1,z2)  =  F(wl5w2)  =  f[01(w1,w2),02(w1,w2)]  , 
by  an  analogous  procedure.  The  Cauchy  inversion  formula  implies  that 


F(W1’W2)=2^I 


f  [h(w1,n),n]{31g2[h(w1,n),n]92h(w1>ri)  +  32g2[h(w1,n)  ,n]  }dn 


g2[h(w1,n) ,n]  -  w2 


(5.4.12) 


However , 

S^g, [h(w  ,n),n]  •  3  h(w  ,n)  =  1 

111  11  (5.4.13) 

31g1[h(w1,n) ,n]  •  32h(w1,n) +  32g1[h(w1,n) ,n]  =  0  . 

Substituting  Equations  (5.4.13)  into  Equation  (5.4.12), 


F(wrw2)  =  2^1 


f  (h,n){-3;Lg23281  +  31g132g2)dn 
(3^)  (g2  ~  «2) 
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(5.4.14) 


where 


h  =  h(w1,n) 

81  =  Si » n)  , n] 

%2  =  §2^^w1,T1^  ' 


But 


-3^2 [h(wi>  n) ,  32si t h(wi»r')  >n]  +  3181[h(w1>n)  ,n]  32g2th(w1,n)  ,n] 


■  J[h(w1,n),n]  , 

i.e.,  the  Jacobian.  Therefore, 

1  f  [h(wx,n)  .nlJthCv/^^.n)  ,n]dn 
f(w1«w2)  -  27ri  a1g1[h(w1,n),n]  (g2[h(w1,n) ,n]-w2>  * 
c’ 


(5.4.15) 


Now  applying  the  Cauchy  inversion  formula  to  the  integrand  in  Equation 
(5.4.15),  i.e.. 


f [h(w  ,n) ,n]J[h(w1,n) ,n] 
31g1[h(w1,n),n](g2[h(w1,n),n]-w2) 


determines 


1 

2tt  i 


3(w  ,w  ) 

dCdn 

[g1(5.n)-w1] [g2(C,n)-w2] 


c 


(5.4.16) 


1 


F(w1#w2)  = 


4t: 


c  c 1 


3(w  ,w  ) 

f(C,n)  3(g<n)  dCdn 

[g1(C,n)-w1] [g2(5»n)-w2] 
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(5.4.17) 


If  f (?,n)  = A(rQ, 6^, 02>t) ,  the  amplitude,  then 


b  =- 

rs  ,  2 
4n 


c  c 


_  3(Y, >Y~) 

A(ro,Yi(C,n),Y2(C,n),T)  —  dgdn 

,  _  V  t+1  .  .  s+1 

Y-j^iC » n)  Y2(5,n) 


(5.4.18) 


Since 


_  3(y1>Y2)  _ 

A(r0,Yi(C,n),  Y2(C,n),t)  ~'3'(g  ny  °  A(ro,C,n,x)  , 


Equation  (5.4.18)  is  equivalent  to  Equation  (5.4.3)  with  c=e,  c’  =  6, 
cf.  [21]. 


CHAPTER  VI 


ASYMPTOTIC  EXPANSIONS  OF  THE  INTEGRALS 


6.1 


Introduction 

Consider  the  integral  determined  in  Chapter  V, 
r  r 


I  = 


^(ro,ei,e2,T>exP{iT  (-ei2  -  6  2^  }dB1d82 


) ) 


(6.1.1) 


When  n=2,  the  stationary  phase  technique  applies  directly.  When 
n^3,  the  classical  technique  must  be  modified  to  determine  the 
asymptotic  expansion.  In  this  chapter,  the  asymptotic  expansion  of 
the  above  integral  is  determined  for  all  n. 


6. 2  Duistermaat1 s  Theorem 

As  was  noted  in  Chapter  I,  there  are  several  approaches  to  the 
classical  stationary  phase  technique.  A  computationally  simple  approach 
derives  from  an  elegant  theorem  of  Duistermaat  [22].  Specifically, 
let  Q(a)  be  a  symmetric,  diagonal  matrix  with  non-zero  eigenvalues, 
and  let  g(x,a,t)  be  a  bounded  function  with  bounded  derivatives  which 
is  analytic  near  the  stationary  points  of  the  phase.  Then,  for  x -*•“>, 
Duistermaat 's  Theorem  determines  the  asymptotic  expansion 

g(x,a,x)exp{ix < Q(a)x,x/2> }dx  ~ 

- 

(~  )n//2  |detQ(a)  |  1/,2exp{i(ir/4)sgnQ(;.) }  E  ~,Dkg(0,a)x  k, 

T  k=0K- 

(6.2.1) 


where 


PWI4J4 1.  WOW?, t 


D  =  4-  <  |detQ(a) 


9x  ’  3x 


sgnQ(a)  denotes  the  number  of  positive  eigenvalues  of  Q(a)  less  the 
number  of  negative  eigenvalues,  and  n  is  the  dimension  of  the  region 
of  integration.  Applying  Equation  (6.2.1)  to  the  integral  in  Equation 
(6.1.1)  determines 


A(ro,8i,62,T)exp{iT(+612+622) IdB-jdBj 
1/2 

~(2)  ^-exp{+i(Tr/4)}  £  f,DkA(7  ,0)T_k 

T  k=(T‘  0 


where 


n  =  i  (+  2  +  29 

2  —  3B2  ’ 


with  the  sign  of  -^r  corresponding  to  the  sign  of  ,  and  it  has  been 
assumed  that  both  31  and  B2  have  the  same  sign  (otherwise,  sgnQ(a)  *  0 
and  exp{ iirsgnQ(a) /4} = 1) .  Duistermaat ' s  Theorem  determines  the  complete 
asymptotic  series  at  any  regular  point. 

6.3  The  Modified  Stationary  Phase  Technique 


When  n^3,  i.e.,  Q(a)  has  one  zero  eigenvalue,  the  theorem  of 
Duistermaat  cannot  be  applied  directly.  The  asymptotic  series  of 
Equation  (6.1.1)  may  yet  be  obtained  by  a  modification  of  Duistermaat ' s 
approach.  Let  Equation  (6.1.1)  be  expressed  as 


^  _  o 

a(ro,B1,B2,T)exp{iT(+B]  +B2n))dB1dB2 


exp{+iT32  }d82  A(r  , B2>T)exp{+iT6,  }d6. 


(6.3.1) 
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The  asymptotic  series  of  the  integral  over  8^  proceeds  from  Equation 

(6.2.1) .  Each  series  coefficient  is  a  function  of  62*  hence,  each  term 
determines  an  integral  over  32-  ^he  asymptotic  series  of  each  integral 
may  be  determined  by  a  modification  of  the  technique  outlined  in 
Chapter  I.  The  complete  asymptotic  series  is  the  sum  of  the  series 

of  the  individual  integrals  over  82* 

More  explicitly,  applying  Duistermaat' s  Theorem,  in  Equation 

(6.2.1) ,  to  the  integral  over  8^  in  Equation  (6.3.1)  obtains 

* 

A(ro,61,  62»'r)exp{+iT812}d81 

1/2  l  i  k_  -  t 

-(7)  exp{+i(ir/ 4) }  I  (7, )  (7)  A,.(r  ,0,8,)t  K  ,  (6.3.2) 

T  k=0  Z  0  1 


where 


A2k(ro’°,62)  " 


32kA(ro,0,B2’x:) 
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2k 


Each  A01  (r  ,0, 
2k  o 

•  T-f 


82)  determines  an  integral  of  the  form 
A2k^ro’0,B2^exp^-iTB2n^dB2  ' 


(6.3.3) 


Since  A(r,8^,82,T)  and  all  its  derivatives  are  bounded  (Chapter  I), 
the  integral  in  Equation  (6.3.3)  may  be  expressed  as 


:(T)A2k(ro,0,82) 


A2k^ro,0,B2^eXM-iTB2  }dS2 


(6.3.4) 


Expanding  A2k(ro>0,82)  about  82 = 0  determines 


1 


A2k^ro,0,62^€xp^— 1t62  ^d®2  =  a0  exp{±i^32  }d&2  +  al  exP{±iT02  ^2d^2  + 

GO  '00  - 

>  —00  —00 

•  •  •  +  an_2  exp{+ix62n}e2n  ^d82  +  ~  exp  t+iTB^}^11  1RA2k(ro,0,62)d62, 


(6.3.5) 


where 


1  s  d 


3  _  _  _ 


aj  =  (F}  A2k(ro*0) 


EA2k(V°-B2>  ’  ”®21””(A2k<ro,°’  ®2>  -  A2k<V°>  "  (fr)^T*2 k(V°»  • 

j=0  J  dB2J 

i.e.,  the  remainder  of  the  Taylor  series  less  a  factor  of  n  1. 

The  integrals 


J^(t)=  exp{+iT62n)62'1d62 


(6.3.6) 


are  determined  by  contour  integration  (Appendix  B) .  A  partial 
integration  of 


i  eKp{±i,62")62"-1RA2k(ro,0,62)d62 


obtains 


^  exp{+ixB  n)B9n  1RA  ,  (r  ,0,6  )dB9  = -r-  (77-exp{+iTB9n})RA  ,  (r  ,0,B9)dB 
n  —  2  2  2k  o  Li.  it  d62  —  L  2k  o  2  2 


I 


00 
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+exp{+iT &2n  ^ro  ’ 0  ’ 13  2^ 


+  — 
IT 


exp{+iTB2n}^[RA2k(ro,0,e2)]dB2  = 


—  oo  —  oo 


exP{+iTB2n}^[RA2k^o.0,B2)]dB2, 


(6. 


where  the  last  equality  follows  from  the  boundedness  of  A(r,p,x)  and 
its  derivatives  [28].  Let  S  be  an  operator  defined  by 


SA2k(ro,0,e2^  *  dB2  ^RA2k(ro’°,62^  ' 


(6. 


Then,  Equation  (6.3.5)  may  be  expressed  as 


I(T)A2k(ro,0,82)  aoJo  +  VJl+  +“n-2Jn-2  irI(T)SA2k(ro’0’ S2) 

(6. 

or,  alternatively,  as  an  operator  equation 


I(r)(i+-rS)=a  J  (t)+;  J  (T)+...+Sn  _(t)  , 

X I  O  o  11  tl—  4  n-  4 


(6. 


where  1  is  the  identity  operator  and  the  ot^  are  the  operators  which 

carry  functions  to  constants.  Because  Equation  (6.3.10)  is  a  power 

*  1 

series,  the  operator  (1+^-S)  has  a  (right)  inverse;  namely 

OO 

E  (+t)~V  .  Thus, 

«.=0 


I(t)  =  (  E  a .J .  (t))(  E  (+ix)"V) 
k=0  k=0 


(6. 


3.7) 


3.8) 


3.9) 


3.10) 


3.11) 


Therefore,  Equation  (6.3.5)  becomes 
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I(T)A2k(ro’°’e2) 


A2k(ro»°»e2*exp*-iT62  *dP2  = 


J At)  Z  (+ix)“VA  (r  ,0,0  )+J  (t)  Z  (+ir)"  Va„.  (r  ,0,0,)  +  .. 
°  l~0  1=0  * 


n-2 


'  ,  .  X  ^  .  /■ 


+  J  (x)  E  (+ii)  S  A_(r  ,0,8  )-(  E  J,(t))[  E  (+ix)  *S*A9.(r 

n_2  1=0  2k  0  2  j=o  3  *=o  2k  ° 


AJj<T>^+lTrt<;ib>Vv°-s2>>  • 


(6 


Thus,  combining  Equations  (6.3.2)  and  (6.3.12)  determines 


w  _ 
A(ro,01,02,T)exp{iT(+01  +02  )}d01d02' 


n-2 


(^)1/2exp{i(TT/4)}  Z  -^r(f)kT  k(  E  J.  (t)[  z  (+iT)  £(-^-R)  V  (r  ,0)]} 
T  k=0K-  1  j=0  3  a=0  dP,2  2k  ° 


(6, 


For  n^3,  Equation  (6.3.13)  is  the  complete  asymptotic  series. 


,0, b2) ] 

.3.12) 


3.13) 


CHAPTER  VII 


CONCLUDING  REMARKS 


7.1  Introduction 

In  this  chapter,  the  complete  algorithm  is  summarized.  An  example 
with  an  explicit  emitter  geometry,  boundary  conditions  on  the  momentum 
and  a  specific  medium  profile  is  presented.  Finally,  some  promising 
extensions  of  this  investigation  are  outlined. 


7.2  Summary  of  the  Algorithm 

We  consider  wave  propagation  in  a  medium  characterized  by  the 
reduced  Helmholtz  wave  equation 


Ailzli  +  jitai  +  Tf(xH(x,y)  -0  .  (1.3.1a) 

3x  3y 

where  x  is  the  depth,  y  is  the  range,  and  f(x)  is  the  profile.  We 
assume  that 


where 


<P(x>y) 


A(x,y,px,py,x)exp{ix$(x,y,px,py)}dpxdpy  , 


A(x,y,p 


,Py,t)  ~ 


Z  \(x, 
k  k 


y»p„ 


(py  )  T 


-k 


(1.3.2a) 

/ 


and  4>(x,y,p  ,p  )  has  the  form 

x  y 

<Kx,y,p  ,p  )  =  xp  +  yp  -  S(p  ,p  )  .  (3.3.7) 

x  y  a  x  x  y 

Let  y  =  g(x)  specify  the  geometry  of  the  emitter  centered  at  (x,y),  with 

the  boundary  momentum,  p =  (p  ,p  ),  and  the  spatial  derivatives  (in  the 

x  y 

x-direction)  of  the  momentum  at  the  emitter  being  known. 
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The  algorithm  begins  by  forming  the  Hamiltonian 
H  =  p  2  +  p  2  -  f  (x)  =  0  . 

*  y 


(3.3.1) 


Then  the  Lagrange  manifold  is  determined  from  Equations  (3.3.2),  (3.3.3), 
and  (3.3.4) 


3S  .-1.  2  2. 

x  =  —  =  f  (p  +p  ) 

3p  x  y 

x 


(7.2.1a) 


a2s  .  as  .  n/  ,  as 

3px3py  x  ~  3py  Py)  ~  3py  * 


(7.2.1b) 


where  O(p^)  is  an  arbitrary  function  which,  for  simplicity,  we  assume 
a  polynomial  in  p^  with  constant  coefficients.  These  coefficients  are 
determined  by  noting  that  on  the  emitter.  Equations  (7.2.1a)  and 
(7.2.1b)  become 


-  3S(p)  2  —  2. 

x  =  r  =  f  (p  +  p  ) 
3px  *x  *y 


(7.2.2a) 


y  y 


(7.2.2b) 


Solving  Equation  (7.2.2a)  for  p^  and  substituting  into  Equation  (7.2.2b) 
determines  an  equation  in  x  and  p^,  e.g.. 


g(x)=G(x,py)  , 


(7.2.3) 


where  x  and.  p^  are  both  known  on  the  emitter.  Successive  differentiations 
of  Equation  (7.2.3)  with  respect  to  x,  evaluated  at  the  emitter,  lead  to 


a  system  of  linear  algebraic  equations  for  the  coefficients  in  0(py)  in 
terms  of  x,  py  and  the  spatial  variations  of  the  momentum  on  the 
emitter.  Solving  for  the  coefficients  determines  S(p  ,p  )  and,  hence, 

x  y 

4>(x,y,p  ,p  ) . 

x  y 
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The  equation  of  the  caustic  in  momentum  space  is  determined  from 

the  Hessian  of  ij>(x,y,p  ,p  ),  i.e., 

x  y 


=  0  ,  (7.2.4) 


and  is  projected  onto  coordinate  space  through  the  Lagrange  manifold. 
Equations  (7.2.1a)  and  (7.2.1b).  That  is,  any  (Px»Py)  satisfying 
Equation  (7.2.4)  determines  a  corresponding  point  (x,y)  in  coordinate 
space  defined  by 


det 


9p  3p 
y  x 


3\ 

S2* 

3P  2 

X 

3px3py 

2 

3d> 

3  (j> 

»P, 


as(p)„  --as(p) 

x  3p  *  y  3p  • 
rx  *y 

The  locus  of  these  points,  (x,y),  defines  the  caustic  in  coordinate 
space. 

The  asymptotic  analysis  of  the  field,  i.e.,  the  integral  in 

Equation  (1.3.2a),  at  any  point  (r0>P0)  proceeds  from  an  analysis  of 

d> (r , p )  at  (r0»P0)>  (Appendix  A).  The  determination  of  the  asymptotic 

Series  of  the  integral  at  (r  ,p  )  requires  that  ()>(r,p)  at  (r  ,p  )  be 

o  o  o  o 

transformed  to  its  canonical  form  (Chapter  V).  The  asymptotic  series  of 
the  transformed  integral 


I 


'  r 

A(ro,61»e2>T)exp{  iT(+812+B2n)  » 


(7.2.5) 


where 

A(7o,6r82»T)  ~  l  VVei’02)x”k  ’ 

with  the  Jacobian  included  implicitly,  requires  the  determination  of 
the  asymptotic  series  for  each  of 
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qpmipipanuii.'iu.)  i  *■  ji'h***'  . . . . |i. 

IWtSVWCVWM, 


V 


ro,g1,e2)exp{iT(+612+B2n) }ciB  de2- 


(7.2.6) 


The  asymptotic  series  of  each  integral  is  determined  in  Chapter  VI.  The 
A^(rQ, 3^,62)  are  determined  from  the  transport  equation  (Chapter  IV). 

The  complete  asymptotic  series  of  integral  in  Equation  (7.2.5)  is  the 
sum  of  the  asymptotic  series  of  each  integral  in  Equation  (7.2.6)  for 
all  k>  0. 


7.3  Example  of  the  Algorithm 

Consider  a  medium  represented  by  a  linear  profile,  f(x)  =x.  We 
investigate  the  far  field,  i.e.,  the  distances  involved  are  much  larger 
than  the  dimension  of  the  emitter,  taken  here  to  be  a  line  source 
centered  at  (5,-4)  with  slope  y = 2x  producing  a  wave  of  constant  amplitude, 
A(r,p)  =1.  Let  the  magnitude  of  the  momentum  at  the  emitter  be  |p|  =5; 
for  definiteness,  let  the  components  of  the  momentum  and  the  spatial 
variation  of  the  momentum  at  (5,-4)  be 


where  the 
the  field 


Pxo 

=  -2.0 

pyo 

=  1.0 

Pxo 

=  -0.67 

py° 

=  -0.83 

Pxo 

=  2.05 

p;’o 

=  2.96 

c 

=  10.2 

pjs 

=  31.9 

(7.3.1) 


primes  indicate  derivatives  with  respect  to  x.  We  represent 
by 


iKx.y)  = 


A(x,y,p  ,p  ,T)exp{iT<j>(x,y,p  ,p  )  }dp  dp  ,  (1.3.2a) 

x  y  x  y  x  y 


A(x,y,px,py,-r) 


I 

k 


Ak(x,y,Px’ 


t 


where 
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i.e.,  integral  in  Equation  (1.3.2a)  is  an  asymptotic  solution  of  the 
reduced  Helmholtz  wave  equation  [Equation  (1.3.1a)]. 

In  this  case,  the  Hamiltonian  [Equation  (3.3.1)]  becomes 

“W-*-0  <7-3-2> 

which  leads  to  the  Lagrange  manifold  [Equations  (3.3.3)  and  (3.3.4)], 


x 


(7.3.3a) 


y  ■t^-  =  2p  p  +  0(p  ) , 

3py  *  y  y 


(7.3.3b) 


For  definiteness,  let  0(py)  be  a  polynomial  in  p^  to  the  third  order. 


i.e. , 


e(Val  +  a2Py  +  a3py2  +  Vy3 


(7.3.4) 


The  a^in  Equation  (7.3.4)  may  be  determined  from  Equations  (7.3.1), 
(7.3.3a,  b)  and  the  emitter  geometry,  y  =  2x.  The  substitution  of  p 


and  y = 2x  into  Equation  (7.3.3b)  leads  to 
2x  =  2p 


/  2  .  ,  2^  3 

/ x  -  p  +  a  +  a_p  +  a  p  +  a.p 
”  1  2  y  3  y  4  y 


(7.3.5) 


y  -y 

Three  differentiations  with  respect  to  x  then  obtain  a  system  of  linear 

equations  which,  with  the  momentum  conditions  at  (5,-4)  specified  by 

Equation  (7.3.1),  determine  a^  =  3,  a2  =  -5,  a^l,  a^  =  l.  Thus,  from 

Equations  (7.3.3)  and  (7.3.4), 

2  2 
x  =  p  +  p 
x  y 


y  =  2pxpy  +  3  -  5Py  +  Py2  +  Py3  , 


(7.3.6) 


and 


4>(x,y,Px,Py)  =xpx  +  yPy 


2  _  5  J 

PxPy  -3py  +  2Py 
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(7.3.7) 


The  equation  of  the  caustic  is  determined  by  equating  the  Hessian 


determinant  of  <ji(x,y,Px»Py)  to  zero,  i.e., 


det 


2P. 


L2py 


2p, 


2p  -  5  +  2p  +3p 
x  y  y  j 


=  0  , 


yielding 

4p  2  -  4p  2  -  lOp  +4p  p  +  6p  p  2  =  0  .  (7.3.8) 

x  y  x  x  y  x  y 

Equation  (7.3.8)  is  the  equation  of  the  caustic  in  momentum  space.  Any 

p =  (p  ,p  )  satisfying  Equation  (7.3.8)  determines  a  corresponding  point 

x  y 

(x,50  in  coordinate  space  defined  by  the  Lagrange  manifold.  Equation  (7.3.6), 


y=2pp  +  3  -  5p  +p2  +  p3  . 

J  x  y  y  y  y 

The  locus  of  these  (x,y)  is  the  equation  of  the  caustic  in  coordinate 
space.  Figure  1. 

The  asymptotic  expansion  of 


i|»(x,y)  = 


A(x,y,px,py,T)exp{iT<{>(x,y,px,py)}dpxdpy 


(1.3.2a) 


proceeds  by  transforming  the  phase  to  its  canonical  form  (Chapter  V) . 
We  illustrate  the  procedure  by  determining  the  first  two  terms  in  the 
asymptotic  series  at  (x,y)  =  (2,2),  which,  following  Appendix  A,  is  the 
cusp  point.  Hence,  at  (2,2)  the  canonical  form  of  4>(x,y,px,py)  is 

M2,2,61,82)  =^>(2,2,px,py)  +  812  + 824=  <f>(2,2,px,py). 
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HORIZONTAL  RANGE  (y)-IN  DIMENSIONLESS  UNITS 


FIGURE  1.  THE  CAUSTIC  CURVE 


I*- 


U,W  i..UUP'-fji  M&IWJWP 
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From  the  Lagrange  manifold  conditions.  Equations  (7.3.6),  the 

corresponding  momenta  at  (x,y)  =  (2,2)  are  (px*py)  =  (1,1).  Expanding 

<Kx,y,p  >P__)  In  a  Taylor  series  at  (2, 2, 1,1)  yields 
x  y 

<H2,2,px,py)  =  y§-  (px-l)2-2(px-l)(py-D-  (Py-1)2 

-i(px-l)3-(pi-l)(py-1)2-f(py-1)3-i(py-l)4  . 

Translating  the  origin  to  (p  ,p  )  =  (1,1)  and  applying  the 

x  y 

principal  axis  transformation,  Equation  (5.2.2),  obtains 


4)(2,2,a1,a2)  12  '  a2  (1  * ^l  +  3a2  +  32^1  ~16ala2  +  64a2  ^ 

,  ,1  Is  14 

2a2(32al _ 4)al  64°!  * 


Then,  the  coordinate  transformation 


Q-/i  5.  1.3  2  1  .  1  2. 

B1  a2(1~8al  3a2  32^1  “  16alCt2  + 64^2  ^ 


1/2 


-1/2 


2,1  _lw.  5  .  1  .  3  2  i  1  2.  .  Q. 

1  (32?1  4)(1~  8al+3a2  +  32al  ~  16ala2  +  64a2  )  (7.3.9) 

_  r,l  1^2m  5.1.321  ,  1  2. _1  1  .1/A 

B2  1^37*1  4)  (1  "  ^l4-  3a2+l2al  l6ala2+64a2  )  -  64} 


carries  <t>(r,p)  at  (2, 2, 1,1)  to 


~  19  2  4 

4(2,2, erS2)  +Bj 


Under  the  coordinate  transformation  in  Equation  (7.3.9),  the 
integral  in  Equation  (1.3.2a)  becomes 


4(2,2) 


A(2,2,61,e2,T)exp{iT(Y|-  612  +  324))dB1dB2  ,  (7.3.10) 


where 
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A(2,2,B1,32,t)  -  E  Ak(2,2,31,62)T-  , 

k 

3(a1>a2) 

with  the  Jacobian  -r-rz — — r-  ,  included  implicitly,  e.g.,  Chapter  V. 

The  determination  of  the  asymptotic  series  proceeds  from  Equations 
(6.3.2)  and  (6.3.5)  and  Appendix  B,  with  n=4.  The  first  two  terms  of 
the  asymptotic  series  at  (x,y)  =  (2,2)  are  then 


<|/(2,2) — 0.77exp{-^p}r(-|)T  ^^{cos  +  isin-^-} 

-0. 54exp{^Li}r (|)  T'5/4{cos  y+  isin  ^  } . 

7 . 4  Suggestions  for  Future  Research 

In  this  investigation,  a  partial  differential  equation  of  classical 
physics,  the  reduced  Helmholtz  wave  equation,  has  been  studied  using 
theorems  from  topology,  more  specifically,  singularity  theory,  in  an 
essential  way.  Singularity  theory  is  a  well-developed  research  area  in 
abstract  mathematics,  e.g.,  [25];  since  the  announcement  of  Thom's 
Classification  Theorem,  singularity  theory  has  been  the  subject  of 
intense  interest  among  pure  and  applied  mathematicians  and  scientists 
[52].  In  the  course  of  this  investigation,  a  number  of  promising 
extensions  of  this  research  have  become  apparent. 

The  mathematics  of  this  investigation  focused  on  the  reduced 
Helmholtz  wave  equation  in  two  spatial  variables.  The  analysis  led  to 
caustics  which  could  be  characterized  by  the  A-series  Thom  potentials. 


An  obvious  generalization  is  to  consider  the  reduced  Helmholtz  wave 
equation  in  three  variables,  which  would  lead  to  caustics  corresponding 
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to  the  D-series  (umbilic)  Thom  potentials.  The  determination  of  the 
caustic  geometry  in  coordinate  space  would  proceed  as  in  Chapter  III 
except  that  the  phase,  (j>(r,p)  would  require  a  Hessian  of  rank  3.  The 

field  would  be  represented  by  an  integral  of  the  form 
fff 


■Kx.y.z)  = 


A(x,y,z,p  ,p  ,p  )exp{iT4>(x,y,z,p  ,p  ,p  )  }dp  dp  dp  , 

a  y  £4  a  y  £  x  y  z 


(7.4.1) 


where 


-k 

A(x,y,z,p  , p  , p  t)~  I  A,  (x ,  y ,  z ,  p  ,  p  ,  p  )  t  . 

Ayz  Ayz 

At  those  field  points  where  the  Hessian  of  (J>(r,p)  was  completely 
degenerate,  the  integral  in  Equation  (7.4.1)  must  be  transformed  to  an 
integral  whose  phase  was  the  appropriate  umbilic  (Appendix  A) .  The 
stationary  phase  technique  developed  above  is  not  valid  at  such  field 
points.  Consequently,  a  further  modification  of  the  technique  would  be 
required  to  determine  the  asymptotic  series. 

The  general  stationary  phase  technique  for  multiple  integrals 
could  itself  be  studied.  As  was  noted  in  Chapter  I,  integral  in 

Equation  (1.2.9)  may  be  transformed  to  an  integral  of  the  form 

' 

exp{ ixu)dg(u)  , 

where  g(u)  is  the  integral  of  A(r,p,t)  over  the  region  <Hr,p)<_u. 

The  asymptotic  behavior  of  the  original  integral  is  now  determined  by 
an  analysis  of  the  singular  points  of  g(u).  It  would  be  interesting 
to  investigate  the  utility  of  Thom's  Classification  Theorem  in  the 
analysis  of  such  integrals. 


66- 

The  physical  phenomenology  considered  in  this  investigation  was 
concerned  with  the  propagation  of  unattenuated  waves  through  a  non- 
dispersive  medium.  An  extension  of  the  above  algorithm  to  include 
attenuated  waves  propagating  in  a  dispersive  medium  seems  possible. 
Indeed,  the  extension  of  the  algorithm  to  include  other  linear  "wave 
equations",  e.g.,  the  Klein-Gordon  equation,  appears  very  promising. 

As  this  investigation  was  concerned  with  classical  physics, 
"semi-classical"  wave  phenomena,  e.g.,  barrier  penetration,  were  not 
considered.  It  appears,  however,  that  the  physical  phenomenology 
motivating  Maslov  and  Arnold  was  primarily  semi-classical  [2,4].  It 
would  be  interesting  to  determine  if  the  above  algorithm  could  be 
extended  to  include  semi-classical  wave  mechanics,  generalizing  the 
results  of  Eckmann  and  Seneor  [24]. 
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APPENDIX  A 


THOM'S  THEOREM  AND  THE  CLASSIFICATION  ALGORITHM 

A.l  Introduction 

Thom's  Classification  Theorem  states  that  a  class  of  functions  <p(x,y), 
noted  below,  may  be  transformed  into  certain  canonical  forms,  commonly 
referred  to  as  the  "Thom  potentials".  The  equivalency  between  <|>(x,y)  and 
its  corresponding  Thom  potential  derives  from  topological  considerations; 
specifically,  one  is  obtained  from  the  other  through  a  diff eomorphism,  a 
differentiable  coordinate  transformation  with  a  differentiable  inverse. 
4>(x,y) ,  then,  is  merely  in  a  different  algebraic  form  than  its  corresponding 
Thom  potential.  In  this  appendix  a  formal  statement  of  Thom's  Theorem 
is  given  and  an  algorithm  that  determines  the  canonical  form  corresponding 
to  a  given  function  is  detailed. 

A. 2  Thom's  Classification  Theorem 

Up  to  the  addition  of  a  non-degenerate  quadratic  form  in  other 
variables  and  up  to  multiplication  by  +1,  an  r-parameter  family  of 


smooth 

functions,  l^_r^_4,  is 

equivalent  up  to  a  diff eomorphism 

to  one 

of  the  following: 

r 

Degenerate  Term 

Universal  Unfolding 

Name 

1 

3 

X 

3^. 

x  +  k^x 

Fold 

4 

4  2  , 

2 

X 

x  +  k^x  +  k£X 

Cusp 

3 

5 

X 

5  3  2 

x  +  kjX  +  k£X  +  k^x 

Swallowtail 

3 

3  2 

x  -  xy 

3  2  2 

x  -  xy  +  k.y  +k„x  +  k.y 

Elliptic 

1  o 

X  4  J 

Umbilic 

3 

J  .  J 

x  +  y 

x  +  y J  +  k.  xy  +  k„x  +  k.y 

X  4  J 

Hyberb  >lic 

6 

6  ,  4  ,  3  ,  2  . 

Umbilic 

4 

X 

x  +  k,x  +  k-x  +  k_x  +  k.x 
12  3  4 

Butterfly 

4 

2,4 
x  y  +  y 

2  4,.  2  2 

x  y  +  y  +  k.x  +  k„y  +k,x  +  k,y 

Parabolic 

L  4  J  H 

Umbilic 

Thom's  Theorem  has  been  rigorously  proven  in  a  series  of  papers  by 
Mather  [36-41]. 
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mined  from 


i.e. , 


32<Kr  )  32$(r  ) 

o  ,  o  _  n 

e  +  — — — -  e„  -  0 


3x 


2  1  3x3y  2 


32<Kr  )  9 2<t>  (r  ) 

- 2_  e  +  - 2_  e  =  0 

3x9y  el+  3y2  e2 


If  =  1,  then 


32(Kr  ) 

o 


3x  e, 

e2~  ~  2  —  ^ 

32*(rQ) 


3x3y 


Then  form 
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f (t)  =  <J>(X0  +  te^  yQ  +  te2). 

The  Thom  potential  corresponding  to  <Kx,y)  at  rQ  is  determined  from  the 
Taylor  series  of  f(t)  at  t = 0;  specifically,  the  first  non-vanishing  term 
in  the  Taylor  series  determines  the  equivalent  Thom  potential.  In  a 
neighborhood  of  rQ  there  exists  a  coordinate  transformation  that  carries 
(|>(x,y)  to  the  degenerate  term  of  its  corresponding  Thom  potential. 

(i)  If  f J(0)  = fIJ(0)  =  0,  f 111 (0)  +  0,  the  corresponding  Thom  potential 
is  the  fold. 

(ii)  If  f*(0)  =  fXI(0)  =  fm(0)  =  0,  fIV(0)  t  0,  the  corresponding  Thom 
potential  is  the  cusp. 

(Hi)  If  f *(0)  =  fI^(0)  =  f***(0)  =  f *V(0)  =  0,  fV(0)  /  0,  the  corresponding 
Thom  potential  is  the  swallowtail. 

(iv)  If  fX(0)  =  f ri(0)  =  fm(0)  =  fIV(0)  =  fV(0)  =  0,  fVI(0)/0,  the 
corresponding  Thom  potential  is  the  butterfly. 

(v)  If  fI(0)  -  fI1(0)  =  fni(0)  =  fIV(0)  -  fV(0)  -  fVI(0)  -  0,  there  is 


no  corresponding  Thom  potential. 
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(It  should  be  noted  the  case  f*(0)  =0,  f^CO)  #0  cannot  occur  here 
2  — 

because  det(d  $(r  ))^0.) 

o 

CASE  III:  d2(j>(ro)=0,  i.e.,  each  3ij<t>(ro>  =  0,  i,j  =  1,2. 

(i)  If  not  only  the  second  derivatives  of  <j>(x,y)  at  rQ  equal  zero, 
but  also  the  third  derivatives,  i.e.,  the  third  order  terms  in 
the  Taylor  series  at  rQ,  there  is  no  corresponding  Thom  potential. 

(ii)  If  the  third  derivatives  of  <j>(x,y)  at  r^  are  not  all  equal  to 

zero,  the  third  order  Taylor  series  may  be  expressed  as 
3  2  2  3 

alx  +a2X  y  +  a3xy  +a4y  »  (A. 3.1) 

where  the  a^  are  the  appropriate  Taylor  series  coefficients 

_  3 

evaluated  at  rQ.  If  a^ ^  0,  dividing  Equation  (A. 3.1)  by  y  and 

equating  to  zero  determines 
3  2 

F(t)=a^t  +  ^2t  +a3t+a4  =  0  * 

where  t  =  x/y.  (If  a.^  =  0  and  a^  f  0,  interchanging  x  and  y 
yields  an  analogous  cubic.)  Four  root  combinations  lead  to 
corresponding  Thom  potentials.  That  is,  for  four  root  combin¬ 
ations  there  exists  a  coordinate  transformation  in  a  neighborhood 


of  r 

}  that  carries  <KO  at  to  the  degenerate  term  on  the 

Thom 

potential: 

(a) 

three  real  equal  roots  : 

fold,  x^ 

(b) 

three  real  unequal  roots  : 

3  2 

elliptic  umbilic,  x  -  xy 

(c) 

three  real  roots,  two  equal  : 

2 

parabolic  umbilic,  x  y  +  y 

(d) 

one  real  root  and  one  complex  : 
conjugate  pair 

3 

hyperbolic  umbilic,  x  +y 

If  a^  *  a^  >*  0,  a2>&2^0  the  corresponding  Thom  potential  is  the 


parabolic  umbilic.  If  a^  =  a^  =  0  and  one  of  or  a^  equal  zero,  there 
is  no  corresponding  Thom  potential. 


APPENDIX  B 


DETERMINATION  OF  THE  CONTOUR  INTEGRALS 


B.l  Introduction 


In  Chapter  VI  the  determination  of  the  asymptotic  series  required  the 
evaluation  of  integrals  of  the  form 


I  = 


exp{iTxn}xmdx  , 


(B.1.1) 


where  m  and  n  are  integers,  n>m  +  l.  These  integrals  are  evaluated 
here. 


B.2  Transformation  to  Common  Form 


Consider  the  integral  in  Equation  (B.1.1).  If  m  is  odd  and  n  is 
even,  then,  by  symmetry 


exp{iTxn}xmdx  =  0  . 


If  m  and  n  are  both  even. 


exp{ iTxn}xmdx  =  2 


exp{iTxn}xmdx  . 


(B.2.1) 


) 


Then  by  a  change  in  variable  to  t = tx  ,  Equation  (B.2.1)  becomes 

sintdt 


exp  { irx11  }xmdx  =  — ^ 

exEiitl  dt=  2  , 

1-k  ac  k1 

costdt  . 
1-k  1 

nr 

t  nx  J 

t 

1-k 


},  (B.2. 2) 


0  0  0 
where  ks  (m+l)/n. 

If  n  is  odd,  express  Equation  (B.2.1)  as 


0 


* 

exp{ iTxn}xmdx  ■ 

exp{  Itx11  }xmdx  + 

J 

co  »a 

°  0 

exp{ iTxn)xmdx. 


(B.2. 3) 


<>  i»  i  n*. .i 


i A. 


j  - . — -*■ 
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Then,  by  a  change  of  variables  to  t  =  -tx  , 
0 


* 

/  1  \  ® 

exp{irx  }x  dx =  ^  { 

' 

costdt 

1-k 

nT 

t 

00  0  C 

sintdt 

1-k 


}. 


Similarly,  by  a  change  of  variables  to  t  =  tx 

O 

sintdt 


1-k 


exp{  ixxn  }xmdx  = 

' 

costdt  . 
1-k  +  1 

J  nc  j 

t  J 

}  . 


0 


B. 3  Evaluation  of  the  Integrals 

The  trigonometric  integrals  in  Equations  (B.2.2),  (B, 
(B.2.5)  can  be  evaluated  by  considering  the  integral 


[  exp{ iz}dz 

1  I  1-k 


J  2 

c 

where  c  is  the  contour  in  Figure  2. 

From  Cauchy's  Integral  Theorem  [35] 


[exp{iz}dz 

rR 

exp{ix}dx  , 

exp{ iz}dz 

J 

1-k 

X 

. 

1-k 

z 

c 

+ 

p  p 

exp{-y}idy 

CR 

exp{iz}dz  _ 

(iy) 


1-k 


1-k 


Consider  first 


exp{ iz}dz 
s 


-R 


(B.2.4) 


(B.2.5) 


2.4)  and 


(B.3.1) 


(B.3.2) 


where  a  *  1  -  k.  Since  z  =  Rexp{ iO} ,  and  hence  dz =  iRexp(io}dO, 
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FIGURE  2.  THE  CONTOUR  OF  INTEGRATION 
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exp{ izldz  _  exp{ iRexp(iQ) }iRexp{ i0}dO 
zS  RSexp{is0} 


=  iR  exp{iR(cos0+isin0) } [cosk0+isink0)d0 


_<  |  iR"'"  S  exp{iR(cos0+isin0)  }  [cosk0+isinkO]d0 )  . 


|  iR^  S  exp{ iR(cosQ+isin0) 1 [cosk0+isink0]d0| <R^  S  ) exp{ -RsinO}d0 


<  R  exp{-2R0/ir}d0  , 


where  the  last  inequality  follows  from  O£0^tt/2,  hence  sin0<iT/2. 

Tf/2  , . 

m+l-n 

R1  S  exp{-2R0/ir}d0  =  ^  R  S (l-exp(-R) )  =  —  R  n  (l-exp{-R}). 


Now  taking  the  limit  as  R  -*■  <*> 
m+l-n 

lim  -^R  n  (  l-exp{ -R})  ■>  0  , 

R^oo 


jxp{  iz } 


dz  =  0  . 


Next  consider 


oxpf iz}dz 


=  r(k). 


i.e.,  the  Gamma  Function,  [1];  hence,  taking  the  limit  of  Equation  (B.3.3) 

as  p+0  and  R -*<*>,  simultaneously,  i.e., 

R  R 

iiml  ^£ilx>dx  =  [cos  (kir/ 2)  +  sinCkir/2)]  SiEkzldZ}  , 

Rx»  X  J  y1-k 

p+0  P  p 


determines 


cqgydx  =  r(k)cos(kTr/2)  =  T(— )cos[(m+l)Tr/2n] 
xl-k  n 

(B.3.4) 


*  r(k)sin(kir/2)  -r(^)sin[(m+l)TT/2nl 
xl  k  n 


Summarizing  from  Equations  (B.2.2),  (B.2.4),  (B.2.5)  and  (B.3.4), 


if  n  is  even 


2r(— ) 

n 


m+l  {cos[  (m+l)7t/2n]  +  isin[  (m+l)ir/2n] }  m  even 

exp{ ixxn}xmdx =  nT  n  , 


m  odd 


and  if  n  is  odd 


2r(— ) 

n 


exp{iTxn}xmdx =  ntn 


2r(S±i) 

n 


cos[ (m+l)w/2n] 


sin[ (m+l) n/2n] 


m  odd 
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